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Abstract 

Using  the  commercially  available  FLUENT  3-D  flow  field  solver,  this  research 
effort  investigated  vortex  breakdown  over  a  delta  wing  at  high  angle  of  attack  (a)  in 
preparation  for  investigation  of  active  control  of  vortex  breakdown  using  steady,  along- 
core  blowing.  A  flat  delta-shaped  half-wing  with  sharp  leading  edge  and  sweep  angle  of 
60°  was  modeled  at  a  —  18°  in  a  wind  tunnel  at  Mach  0.04  and  Reynolds  number  of  3.4  x 
1 05.  A  hybrid  (combination  of  structured  and  unstructured)  numerical  mesh  was 
generated  to  accommodate  blowing  ports  on  the  wing  surface.  Results  for  cases  without 
and  with  along-core  blowing  included  comparison  of  various  turbulence  models  for 
predicting  both  flow  field  physics  and  quantitative  flow  characteristics.  FLUENT 
turbulence  models  included  Spalart-Allmaras  (S-A),  Renormalization  Group  k-s, 
Reynolds  Stress  (RSM),  and  Large  Eddy  Simulation  (LES),  as  well  as  comparison  with 
laminar  and  inviscid  models.  Mesh  independence  was  also  investigated,  and  solutions 
were  compared  with  experimentally  determined  results  and  theoretical  prediction.  These 
research  results  show  that,  excepting  the  LES  model  for  which  the  computational  mesh 
was  insufficiently  refined  and  which  was  not  extensively  investigated,  none  of  the 
turbulence  models  above,  as  implemented  with  the  given  numerical  grid,  generated  a 
solution  which  was  suitably  comparable  to  the  experimental  data.  Much  more  work  is 
required  to  find  a  suitable  combination  of  numerical  grid  and  turbulence  model. 
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COMPUTATIONAL  FLUID  DYNAMICS  INVESTIGATION  OF  VORTEX 


BREAKDOWN  FOR  DELTA  SHAPED  WING  AT  HIGH  ANGLE  OF  ATTACK 

I.  Introduction 

While  physical  experimentation  gives  more  accurate,  realistic  results  which 
include  nonlinear  effects,  it  is  often  laborious,  costly  and  time-consuming,  and  this  is 
where  computational  fluid  dynamics  (CFD)  and  numerical  modeling  and  testing  enter  the 
picture.  While  computer  processors  will  likely  never  be  fast  enough  to  meet  everyone’s 
desires,  today’s  numerical  processing  capability  lends  itself  to  more  complex  solvers  and 
accordingly  more  accurate  solutions,  with  lower  cost  and  eventually  less  labor  and  time. 

This  research  effort  consists  of  using  commercial  software  to  generate  a 
numerical  mesh  and  flow  field  model  which  accurately  simulates  and  provides 
quantitative  and  qualitative  predictions  comparable  to  results  obtained  from  wind  tunnel 
testing  of  a  half  delta  wing  at  high  angle  of  attack  and  low  Mach  number.  This  is  part  of 
a  larger  research  effort  to  possibly  eliminate  the  need  for  conventional  control  surfaces  - 
namely  ailerons,  flaps,  elevator,  and  rudder  -  by  providing  closed-loop  control  of  a  delta 
wing’s  vortex-dominated  flight  dynamics;  this  is  done  specifically  by  using  jets  on  the 
wing  surface  to  control  vortex  breakdown  (discussed  later  in  this  chapter).  CFD  results 
may  contribute  to  development  of  that  closed-loop  control  system  and  substantially 
decrease  efforts  in  wind  tunnel  testing.  Such  control  potentially  leads  to  increased  lift, 
reduced  drag,  attached  flow,  and  more  favorable  pressure  gradients  on  the  wing  surface  at 
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higher  angles  of  attack,  all  of  which  enhance  aircraft  performance  (particularly  for  fighter 
aircraft)  and  extend  aircraft  structural  longevity.  Fewer  control  surfaces  also  result  in 
reduced  size,  weight,  and  radar  cross-section,  among  other  benefits.  Such  benefits  have 
pertinence  and  obvious  value  to  the  warfighter  (Gutmark  et  al,  2000).  This  thesis  does 
not  include  investigation  of  different  types  of  vortex  breakdown  control  but  focuses 
rather  on  the  method  of  choice  for  control  jets  -  along-core,  steady  blowing. 

This  research  effort  works  in  concert  with  current  wind  tunnel  testing  at  the 
University  of  Cincinnati  (UC),  Ohio,  from  which  aerodynamic  data  describing  the  flow 
field  and  evolution  of  vortex  breakdown  on  an  instrumented  delta-shaped  half-wing  test 
article  were  obtained.  Control  jets  along  the  mean  trajectory  of  the  primary  vortex  have 
been  fixed  to  give  system  control  (along-core  blowing).  For  this  testing,  the  60-degree 
swept  wing  was  fixed  at  approximately  15  degrees  angle  of  attack  and  Reynolds  number, 
Re  =  3.4  x  105  (based  on  root  chord  length).  Effect  of  continuous  blowing  into  the  vortex 
core  has  been  characterized  and  quantified  for  different  blowing  momenta.  Periodic 
blowing  is  still  under  investigation. 

Collaborative  parties  include  UC,  Ohio  University  (Athens),  Air  Vehicles 
Directorate  of  the  US  Air  Force  Research  Laboratory  or  AFRL/VACA  (Wright-Patterson 
Air  Force  Base,  Ohio  (WPAFB)),  Dayton  Area  Graduate  Studies  Institute  or  DAGSI 
(Kettering,  Ohio),  and  Air  Force  Institute  of  Technology  or  AFIT  (WPAFB). 

Objectives  and  Scope 

The  three  objectives  of  this  research  effort  were:  first,  develop  a  three- 
dimensional  numerical  mesh  and  flow  model  (specifying  turbulence  model,  order  of 
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accuracy,  unsteady  or  steady  flow  assumption,  and  other  model  parameters)  which 
adequately  and  accurately  represent  the  physical  model  and  wind  tunnel  test  data  and  are 
simple  enough  to  limit  the  amount  of  computation  time  for  obtaining  a  solution;  second, 
generate  numerical  data/solutions  which  correlate  as  much  as  possible  with  the 
experimental  data;  and  third,  vary  parameters  such  as  jet  angle,  jet  location,  jet 
momentum,  wing  angle  of  attack,  and  freestream  velocity,  to  assess  vortex  breakdown 
control  sensitivity  and  optimization.  The  CFD  arm  of  this  project  was  intended  to 
enhance,  not  replace,  the  physical  wind  tunnel  experimentation. 

To  summarize  the  results  of  this  CFD  research  effort,  a  numerical  mesh  was 
developed  and  shown  to  be  adequate  though  not  optimal,  in  order  to  minimize  computer 
processing  time  at  this  stage  of  research.  While  some  mesh  optimization  was  perfonned, 
more  refinement  is  necessary  for  greater  solution  fidelity.  Using  this  sub-optimal 
numerical  mesh  and  carefully  selected  boundary  conditions  and  solver  parameters,  none 
of  the  flow  models  investigated  in  this  study,  excepting  the  Large  Eddy  Simulation  model 
for  which  the  computational  mesh  was  insufficiently  refined  and  which  was  not 
extensively  investigated,  predicted  quantitatively  acceptable  proximity  with  experimental 
data  or  qualitatively  accurate  representation  of  flow  physics.  Numerical  modeling  of 
vortex  breakdown  control  using  along-core  blowing  was  the  primary  objective  of  this 
study.  While  along-core  blowing  was  shown  to  have  an  effect  in  a  two-dimensional  CFD 
model,  attempts  to  achieve  similar  effects  in  the  three-dimensional  case  with  these  grid 
and  turbulence  model  combinations  were  unsuccessful.  Further  investigation  is  required 
to  find  an  optimal  combination  of  grid  and  turbulence  model  for  the  numerical  modeling 
of  steady  along-core  blowing  jets  and  vortex  breakdown  control. 
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Delta  Wing  Background 

Leading-edge  separation  and  vortex  generation  are  aerodynamic  characteristics  of 
flow  over  a  delta-shaped  wing  with  an  angle  of  attack  greater  than  5  degrees,  sharp 
leading  edge,  and  no  camber  (Rusak,  Lamb,  1998:  2).  Boundary  layers  from  the 
windward  (lower)  and  leeward  (upper)  wing  surfaces  separate  at  the  leading  edge,  and 
these  shear  layers  roll  into  a  primary  or  leading-edge  vortex  pair  above  the  upper  wing 
surface,  shown  as  item  1  in  Figure  1.1.  For  a  wing  of  geometry  similar  to  that  used  for 
this  study,  the  boundary  layers  were  observed  to  shed  at  a  frequency  of  5  Hz  (Mitchell, 
Delery,  2001:  409). 

The  leading-edge  vortex,  which  is  characterized  by  high  velocities  and  low  static 
pressure,  increases  in  diameter  and  intensity  as  the  core  follows  a  path  downstream  and 


Figure  1.1  -  Flow  over  a  Highly  Swept  Delta  Wing  with  Sharp  Leading  Edge 


(Hoemer,  Borst,  1985:  18.3) 
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inboard  at  an  angle  slightly  greater  than  the  sweep  angle.  As  the  wing’s  angle  of  attack 
increases,  the  vortex  axial  and  rotational  velocities  increase  and  the  vortex  core  height 
above  the  wing  increases  and  the  core  moves  inboard  (Hoerner,  Borst,  1985:  18.2-3; 
Ekaterinaris,  Schiff,  1990:1).  Primary  vortex  generation  is  nearly  independent  of 
Reynolds  number  due  to  the  miniscule  effective  length,  or  radius  of  curvature,  of  the 
leading  edge;  however,  high  Reynolds  number  flow  does  decrease  vortex  diameter 
because  it  effectively  adds  energy  and  velocity  to  the  core  resulting  in  a  more  tightly 
wrapped  core  (Mitchell,  Delery,  2001:  388;  Novak,  Sarpkaya,  1999:  831). 

The  primary  vortex  pair  creates  lateral,  outboard  boundary-layer  flow  on  the  wing 
surface,  which  collides  with  the  primary  separation  and  results  in  additional  separation 
and  a  corresponding  secondary  vortex  pair.  This  lateral  flow  and  secondary  vortex  pair 
are  shown  as  items  2  and  4  in  Figure  1.1.  The  secondary  vortex  pair  is  weaker,  is  located 
outboard  of  and  rotates  in  a  direction  opposite  to  the  primary  vortices.  Unlike  the 
primary  pair,  the  secondary  vortex  pair’s  strength  and  size  are  dependent  on  Reynolds 
number.  Strength  of  the  secondary  vortex  is  then  a  function  of  area  covered  by  and 
velocity  of  the  lateral  boundary  layer  flow  (Hoerner,  Borst,  1985:  18.2;  Mitchell,  Delery, 
2001:  389). 

Delta  Wing  Lift 

While  slender  delta  wings  with  sharp  leading  edges  have  agreeable  perfonnance 
characteristics  in  supersonic  flight,  the  highly  swept  delta  wing  also  generates  additional 
-  albeit  nonlinear  -  lift  at  high  angles  of  attack  and  subsonic  speeds  due  to  vortex 
generation  and  the  resultant  lower  pressure  on  the  leeward  surface.  For  subsonic  flow,  a 
combination  of  potential  and  vortex  lift  then  comprises  the  delta  wing’s  total  lift. 
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Potential  lift  may  be  detennined  by  applying  linear  lifting  surface  theory, 
assuming  that  flow  remains  attached  over  a  sharp  edge,  as  it  does  over  a  round  edge.  The 
potential  lift  coefficient  is 

CLP  =  Kp -since -cos2  a  (1.1) 


where  a  is  angle  of  attack,  and  Kp  is  constant  of potential  lift  defined  by 


Kp  = 


2-b-T 
S  ■  •  sin  a 


where  b  is  wing  span,  S  is  wing  area,  Too  is  freestream  velocity,  and  /’  is  effective 
circulation  or  ratio  of  rotational  and  axial  velocity  (Kuethe,  Chow,  1998:  498-500). 

As  for  vortex  lift,  Polhamus’  theory  states  that  vortex  formation  causes  the 
stagnation  line  to  move  from  the  wing’s  leading  edge  to  its  leeward  surface,  resulting  in  a 
suction  force,  which  then  increments  the  net  lifting  force  (Polhamus,  1971).  The  vortex 
lift  coefficient  is  defined  by 


.  2  ^  wi  _ _  .  2  cosa 

v  =  Ky  •  sin  a-cosa=  1 - -  Kp  •  sin"  a - (1.3) 

^  F^-sin  a)  cos  A 

where  Ky  is  constant  of  vortex  lift  (and  can  be  detennined  by  the  relationship  in  Equation 
1.3),  Wi  is  induced  velocity,  and  A  is  delta  wing  sweep  angle  (Kuethe,  Chow,  1998:  500). 
Total  lift  coefficient  for  a  delta  wing  becomes  the  sum  of  the  potential  and  vortex  lift 
coefficients. 

Cl  =  CLP  +  CLV  ( 1 .4) 
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Polhamus’  theoretical  curves  hold  only  for  a  freestream  Mach  number,  M&  ~  0. 
The  corrected  Kp  and  Ky  values  for  a  subsonic  Mach  number  become 


(Hoemer,  Borst,  1985:  18.19).  Compressibility  effects  may  be  neglected  for  this  study, 
where  /VC.  =  0.0445  (for  Voo  =  15.4  m/s).  Polhamus’  Kp  and  Ky  curves,  generated  from 
Equations  1.2  and  1.3,  provide  a  reasonably  accurate  analytical  prediction  of  lift 
coefficient,  but  they  overestimate  for  angles  of  attack  greater  than  five  or  six  degrees,  due 
primarily  to  a  phenomenon  known  as  vortex  breakdown,  which  is  addressed  in  the  next 
section  (Guillot,  1999:  8-9).  The  overestimation  also  comes  from  an  assumption  that  the 
vortices  stream  perpendicular  to  the  wing’s  trailing  edge,  which  is  a  valid  approximation 
only  for  delta  wings  with  sweep  angle  greater  than  65  degrees  (Wentz,  Kohlman,  1971: 
159). 

Vortex  Breakdown 


As  first  identified  in  the  1950s,  highly  swept  delta  wings  lose  lift  once  they 
exceed  a  critical  angle  of  attack,  due  to  vortex  breakdown.  As  previously  noted  for 
increasing  angle  of  attack,  vortex  energy  and  velocities  increase  and  vortex  static 
pressure  decreases  (while  total  pressure  increases),  resulting  in  greater  lift;  however,  as 
with  all  good  things,  it  has  its  limitations.  At  the  critical  angle  of  attack,  inner  core 
rotation  has  been  observed  to  reach  about  1,000  Hz,  at  which  point  the  vortex  bursts  or 
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experiences  breakdown  -  rapidly  expanding  in  diameter  and  decelerating  axially  and 
rotationally  (Novak,  Sarpkaya,  1999:  825).  This  breakdown  causes  a  decrease  in  lift  due 
to  increased  static  pressure  and  because  the  vortex  effectively  separates  from  the  wing 
surface  and  flow  reattachment  fails  (Hoemer,  Borst,  1985:  18.15). 


Figure  1.2  —  Vortex  Breakdown  Visualization 
(Office  National  d’Etudes  et  Recherches  Aeronautiques) 


The  flow  field  about  the  delta  wing  now  assumes  the  following  three  divisions. 
First,  approach  flow  is  laminar  and  approximately  irrotational.  Second,  flow  within  the 
breakdown  region  is  typically  stagnated,  reversed,  then  restored  to  its  original  direction 
(for  Reynolds  number  less  than  6,000)  with  large  fluctuations  in  velocity;  it  transitions  to 
turbulent  flow,  and  its  size  has  been  observed  to  be  about  five  vortex-core  diameters  in 
length  and  several  diameters  across.  Third,  a  follow-on  vortex  structure  with  a  larger 
core  radius  continues  downstream.  Increasing  angle  of  attack  beyond  the  critical  angle 
causes  the  breakdown  zone  to  advance  toward  the  wing  apex.  Reynolds  number  has 
virtually  no  effect  on  the  occurrence  of  vortex  breakdown,  but  it  does  affect  the 
breakdown’s  form.  (Rusak,  Lamb,  1998:  2;  Leibovich,  1978:  221-223;  Mitchell,  Delery, 
2001:  386;  Faler,  Leibovich,  1977:  1385) 
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Six  distinct  breakdown  forms  have  been  observed  and  classified  (Faler, 
Leibovich,  1977).  More  recent  research  has  shown  that  five  of  these  forms  are  transient, 
unstable  stages  of  the  most  stable  and  repeatable  spiral  form,  shown  in  Figure  1.3a. 
Further,  vortex  breakdown  is  being  redefined,  or  more  generalized,  as  “the  transformation 
of  a  slender  vortex  into  three-dimensional  forms....  Neither  a  stagnation  point,  nor  a 
region  of  reversed  flow,  nor  the  bridging  of  laminar-turbulent  states  is  necessary;”  neither 
is  the  flow  axisymmetric  or  laminar  (Novak,  Sarpkaya,  1999:  825,  833).  Reversed  flow, 
which  causes  the  bubble-type  breakdown,  shown  in  Figure  1.3b,  and  which  causes  the 
presence  of  two  internal  cells,  depicted  in  the  longitudinal  cross-section  of  a  burst  vortex 


(b) 


Figure  1.3  —  Photographs  of  Vortex  Breakdown  in  Cylindrical  Water  Tunnel: 
(a)  Spiral  Form,  (b)  Bubble  Form  (Leibovich,  1978:  222) 
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in  Figure  1.4,  are  not  observed  in  high-Reynolds-number  cases  (Re  >  3  x  105),  though 
this  structure  was  predicted  by  the  numerical  models  used  in  this  investigation.  A 
stagnation  point  in  the  axial  flow  may  or  may  not  be  present,  and  when  present  it  rotates 
around  the  centerline.  Bubble-type  breakdown  and  four  other  types  all  give  way  to  the 
spiral-type  vortex  breakdown  at  high  Reynolds  number  (Novak,  Sarpkaya,  1999). 


Figure  1.4  -  Time -Average  Streamlines  in  Longitudinal  Slice  of  Bubble  Form  Vortex 
Breakdown,  Revealing  Two-Celled  Structure  (Leibovich,  1978:  230) 

Flow  with  vortex  breakdown  is  inherently  unsteady  -  in  vortex  frequency,  in 
vortex  breakdown  or  stagnation  point  position,  in  vortex  centerline  location,  in  primary 
vortex  winding,  and  in  breakdown  form  or  type.  While  the  boundary-layer  shedding 
frequency  may  remain  relatively  steady,  core  frequencies  vary  in  a  seemingly  random, 
not  periodic,  nature.  Decreased  frequency  is  conducive  to  breakdown,  and  randomness  in 
vortex  frequency  leads  then  to  shifting  of  both  breakdown  and  centerline  locations. 
Variations  in  these  positions  range  from  one  to  several  vortex  core  radii,  and  vortex 
breakdown  position  rotates  about  the  centerline  at  a  frequency  of  about  2  Hz  (for  low 
Reynolds  number).  Primary  vortex  winding  has  been  observed  in  most  cases  to 
arbitrarily  reverse  sense.  For  low  Reynolds  number,  the  spiral  fonn  was  observed  to 
change  to  bubble  form  after  3-5  minutes,  then  back  again  to  a  spiral,  further  indicating 
the  transient  and  unsteady  nature  of  the  bubble  breakdown  type.  Additional  unsteadiness 
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comes  from  rotation  of  the  entire  vortex/vortex  breakdown  fonn,  where  it  rotates  at  about 
1.5  Hz  (observed  for  Re  =  3,120);  that  frequency  changes,  however,  with  changes  in 
rotational  or  swirl  velocity  of  the  vortex,  which  is  also  unsteady.  For  high  Reynolds 
number  (Re  >  2.3  x  105),  the  core  rotates  sufficiently  rapidly  that  the  human  eye  cannot 
discern  the  spiral  breakdown  but  rather  sees  a  conical  shape.  Thus  time-averaged 
solutions  cannot  accurately  depict  every  aspect  of  the  highly  unsteady  physics  of 
instantaneous  vortex  breakdown  flows.  (Novak,  Sarpkaya,  1999;  Leibovich,  1978;  Faler, 
Leibovich,  1977) 

Prediction  of  Vortex  Breakdown  Position.  “The  embarrassing  number  of 
different  theoretical  notions  has  not,  it  must  be  admitted,  led  to  satisfactory  understanding 
of  the  flows  observed”  (Faler,  Leibovich,  1977:  1385).  Research  and  hypothesizing  in 
the  ensuing  25  years  have  not  refuted  this  statement,  though  vortex  breakdown  continues 
to  become  better  characterized. 

Numerous  models  have  been  developed  to  predict  location  of  vortex  breakdown, 
and  this  is  not  a  trivial  task,  given  the  numerous  unsteady  and  nonlinear  conditions 
described  above  (Rusak,  Lamb,  1998:  2;  and  numerous  numerical  studies  described  and 
cited  by  Rusak  and  Lamb).  For  this  study,  the  author  has  chosen  to  use  a  model  proposed 
by  Lance  Traub  that  relies  heavily  on  empirical  data  rather  than  theoretical 
understanding,  since  no  one  has  yet  developed  a  theoretical  model  to  account  for  all,  or 
even  the  most  significant,  unsteady  effects  in  vortex  breakdown. 

Traub ’s  method  uses  data  from  thirteen  different  experimental  studies  with  delta 
wings  of  various  sweep  angles  and  at  various  angles  of  attack  and  provides  equations  to 
curve-fit  the  data.  The  simple  model,  which  consists  of  stepping  through  three  algebraic 
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equations,  predicts  vortex  breakdown  location  within  the  scatter  of  data  for  four  different 
sweep  angles  (Traub,  1996).  The  model  is  based  on  sweep  angles  of  65,  70,  75  and  80 
degrees,  and  thus  does  not  necessarily  provide  accurate  prediction  for  A  =  60°  (as  in  this 
study);  however,  the  model  suffices  for  purposes  of  approximation  considering  that 
precise  identification  of  breakdown  location  remains  largely  subjective  (generally  a 
visual  observation). 

Control  of  Vortex  Breakdown 

Some  key  observations  about  delta  wing  vortices  and  vortex  breakdown,  which 
pertain  particularly  to  control  of  the  phenomenon,  include  the  following.  Increased  angle 
of  attack  and  larger  sweep  angle  produce  greater  rotational  velocity  and  when  the  ratio  of 
swirling  to  axial  velocity  exceeds  about  1.3  breakdown  occurs.  Decreased  rotational 
velocity  (from  decreased  angle  of  attack)  and/or  increased  axial  velocity  of  the  vortex 
cause  breakdown  to  propagate  downstream  of  the  wing  apex,  or  even  delay  breakdown. 
Pressure  reduction  also  stabilizes  the  vortex  (Faler,  Leibovich,  1977:  1398;  Mitchell, 
Delery,  2001:  386).  Thus,  ways  to  inhibit  or  delay  vortex  breakdown  include  decreasing 
angle  of  attack,  increasing  wing  sweep  angle,  or  increasing  swirling  and  axial  velocities 
forward  of  breakdown  position  and  within  the  constraining  ratio  noted  above.  Converses 
of  these  abet  vortex  breakdown. 

Controlling  vortex  breakdown,  whether  delaying  or  encouraging  it,  has  numerous 
positive  effects  on  aircraft  handling  and  performance.  By  delaying  or  preventing  it,  body 
lift  increases,  drag  decreases,  flow  remains  attached  (or  re-attaches),  favorable  pressure 
gradients  abound,  structural  cyclic  fatigue  abates,  and  aircraft  stability  is  augmented.  By 
provoking  vortex  breakdown  on  one  wing,  the  resultant  asymmetry  may  enhance  the 
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maneuvering  performance  of  a  combat  aircraft  (Mitchell,  Delery,  2001:  387-388).  For 
these  reasons,  a  significant  number  of  control  studies  have  been  completed  over  the  past 
four  decades. 

Previous  and  Current  Work  on  Vortex  Breakdown  Control.  Again,  since  this 
effort  does  not  include  investigation  of  different  types  of  breakdown  control,  following  is 
a  summary  of  methods  used,  focusing  more  detail  on  the  method  of  choice  for  control  jets 
-  along-core,  steady  blowing. 

Two  general  categories  of  vortex  breakdown  control  are  via  mechanical  structures 
and  via  pneumatic  methods.  Anthony  Mitchell  and  Jean  Delery  provide  an  excellent 
summary  of  research  efforts  in  their  recent  article  in  Progress  in  Aerospace  Sciences 
(Mitchell,  Delery,  2001).  Mechanical  structures  include  “strakes,  canards,  fillets, 
leading-edge  extensions,  flaps  and  vortex  fences,”  while  pneumatic  methods  include 
steady  and  periodic  “spanwise  blowing,  tangential  blowing,  leading-edge  blowing,  along- 
the -vortex-core  blowing,  trailing-edge  blowing,  leeward  surface  suction,  leading-edge 
suction,  and  suction  along  the  vortex  core”  (390).  While  mechanical  devices  are 
generally  more  robust,  they  add  weight  and  drag;  and  while  pneumatic  techniques  have 
shown  greater  benefit  than  mechanical  devices,  they  are  more  subject  to  contamination 
(particles  clogging  the  injection/suction  ports)  (391-395). 

Suction  and  blowing  both  provide  the  same  result  -  to  reduce  static  pressure  along 
the  vortex  core,  increasing  its  stability  -  though  they  accomplish  it  differently.  Studies 
have  shown  that  each  suction  or  blowing  technique  enhances  aerodynamic  performance 
of  the  wing,  but  “none  of  these  techniques  has  clearly  demonstrated  a  superior  efficiency 
or  effectiveness  in  controlling  either  the  vortical  flow  structure  or  the  vortex  breakdown 
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location”  (415).  Nevertheless,  the  more  promising  of  the  pneumatic  techniques  and  more 
subject  to  current  research  are  along-core  and  periodic  blowing  (395-416). 

Along-core  blowing,  whose  magnitude  is  a  function  of  blowing  mass  flow  rate 
and  freestream  velocity,  adds  momentum  to  the  vortex  core  and  increases  both  axial  and 
rotational  velocities,  allowing  for  a  more  stable  pressure  gradient  along  the  wing  surface 
and  for  a  more  steady  vortex  core.  In  most  investigations,  steady  blowing  has  required  a 
great  deal  of  energy  to  affect  the  vortex  breakdown  location,  but  a  recent  study  at 
Louisiana  State  University  (LSU)  has  shown  that  breakdown  location  manipulation  may 
be  accomplished  with  smaller  flow  rates  (Guillot,  1999).  For  even  less  mass  addition,  as 
with  pulsed  blowing,  vortex  breakdown  delay  and  lift  augmentation  have  been 
demonstrated.  Since  periodic  blowing  appears  to  be  most  effective  at  the  natural 
shedding  frequency  of  the  shear  layers,  difficulty  lies  in  identifying  and  matching  this 
unsteady  characteristic  (Mitchell,  Delery,  2001:  409-411,  415). 

Figure  1.5  shows  flow  visualization  results  from  the  LSU  study  without  (column 
a)  and  with  along-core  steady  blowing  (column  b)  from  a  surface  port  at  x/c  =  0.30  and 
angled  to  intersect  the  vortex  core  centerline  (Guillot,  1999:  54).  In  column  a  of  the 
figure,  the  vortex  core  expands  or  bursts  aft  of  x/c  =  0.30,  and  the  core  progressively 
becomes  more  turbulent  downstream  of  that  point;  column  b  shows  the  effect  of  along- 
core  blowing  into  the  vortex,  where  the  vortex  core  maintains  its  integrity. 
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Figure  1.5  —  Transverse  Cross  Sections  of  Primary  Vortex  for  (a)  No  Blowing  and  (b) 
Along-Core  Steady  Blowing  from  Port  at  x/c  =  0.30  (Guillot,  1999:  54) 

Computational  Fluid  Dynamics  Background 

Conservation  of  mass,  momentum  and  energy  constitute  the  equations  of  fluid 
motion,  and  the  Navier-Stokes  equations  represent  those  conservation  laws  in  partial 
differential  fonn  (Hoffman,  Chiang,  2000.1:  274).  These  partial  differential  equations,  in 
integral  fonn,  are  then  approximated  as  finite-volume  (FV)  expressions  and  refonned  into 
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algebraic  equations  to  allow  for  numerical  computation  within  a  specified 
physical/computational  domain  (Hoffman,  Chiang,  2000.1:  358).  To  reduce 

computational  time  and  complexity,  simplifications  of  these  governing  equations  can 
include  assuming  one-  or  two-dimensional  flow,  inviscid  flow,  steady  flow, 
incompressible  flow,  and/or  first-order  solution  accuracy.  For  this  study,  cases  of  interest 
included  preliminary  use  of  simplified  forms,  but  ultimately  focused  on  three- 
dimensional,  viscous,  second-order  accurate,  compressible,  unsteady  flow  -  as  vortices 
and  vortex  breakdown  are  inherently  unsteady. 

Analysis  of  fluid  flow  using  computational  fluid  dynamics  (CFD)  is  an  iterative 
process  consisting  of  three  basic  steps:  mesh  generation,  numerical  modeling  and 
computation,  and  solution  analysis  and  evaluation. 

Mesh  Generation 

Mesh  or  grid  generation  consists  of  creating  a  set  of  grid  points  along  the 
boundaries  and  throughout  the  domain  of  interest.  For  simple  three-dimensional 
geometries,  a  structured  mesh  may  be  generated,  where  all  volumetric  cells  are 
hexahedrons.  Often  the  physical  domain  involves  geometry  which  is  not  rectangular, 
such  that  coordinate  transformation  must  be  performed  to  convert  to  a  computational 
domain  which  is  rectangular  for  a  structured  mesh  approach.  However,  since  few 
problems  involve  such  simple  geometry,  an  unstructured  mesh  approach  is  often 
preferred,  where  volumetric  cells  may  be  tetrahedrons,  pyramids,  or  any  three- 
dimensional  polyhedrons,  and  where  the  domain  is  already  in  a  computationally  suitable 
form,  such  that  no  transformation  is  necessary  (Hoffman,  Chiang,  2000.11:  356-357).  In 
some  cases  and  for  this  effort,  structured  and  unstructured  regions  of  a  domain  may  be 
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combined  into  a  hybrid  mesh.  While  unstructured  meshes  require  greater  computational 
effort  and  memory,  fewer  cells  are  generally  required  within  a  given  domain. 

Numerical  Modeling  and  Computation 

Once  a  suitable  mesh  has  been  generated,  a  computational  algorithm  numerically 
solves  the  FV  equations  for  fluid  values  at  each  volumetric  cell  center,  from  which  these 
values  can  be  interpolated  to  cell  faces  (FLUENT,  2001:  22.2).  Computations  may  cease 
once  the  solution  has  converged  or  has  come  to  fully  developed  flow.  Normally, 
convergence  is  achieved  when  the  FV  approximation  approaches  the  partial  derivative 
solution,  and  may  be  detected  when  error  residuals  (function  of  difference  between 
previous  and  current  computed  values)  are  sufficiently  reduced,  when  some  integrated 
value  (such  as  lift  coefficient)  becomes  steady  for  a  steady-state  solution,  and/or  when 
flux  (such  as  mass  flow)  is  conserved  (Hoffman,  Chiang,  2000.1:  26;  FLUENT,  2001: 
22.16). 

Solution  accuracy  and  convergence  are  a  function  of  temporal  step  size,  flow 
model  (for  example,  order  of  accuracy,  under-relaxation,  viscous  model,  and  equation 
coupling),  and  mesh.  Steady-state  solutions  converge  more  quickly  with  local  time¬ 
stepping,  where  step  size  (At)  is  determined  and  updated  based  on  local  values  within 
each  FV  of  the  domain.  Local  time-stepping  also  maintains  greater  solution  stability 
within  the  domain.  However,  time  accurate  or  unsteady  solutions  must  use  either  a 
global  time  step  or  a  combination  of  global  and  local  time  steps,  where  the  time  step  must 
lend  stability  to  the  numerical  scheme  and  suit  the  physical  requirements  related  to  the 
problem  (such  as  to  not  exceed  the  vortex  frequencies  if  they  are  to  be  accurately 
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evaluated  in  time).  Time  step  may  be  defined  as  either  one  global  value  or  by  the 
following, 

At 

At  =  (CFL)  ■  — —  (1.7) 

/hnax 

where  Ax  is  spatial  step  size,  Xmax  is  maximum  of  local  eigenvalues  (function  of  flow 
velocities  in  three  dimensions,  fluid  density,  and  speed  of  sound  in  air;  refer  to  FLUENT 
User  Guide  for  further  detail),  and  CFL  is  Courant-Friedrichs-Lewy  condition  (a 
specified  value).  In  the  combined  case  (used  for  unsteady  solutions  in  this  study),  the 
solution  is  driven  to  convergence  at  each  global  time  step  or  physical  time  level. 
(Hoffman,  Chiang,  2000.11:  146-149,  266;  FLUENT,  2001:  22.4.3-4) 

Flow  model  order  of  accuracy  in  space  and  time  also  affects  the  solution 
accuracy.  While  the  governing  equations  of  fluid  mechanics  may  be  discretized  only 
spatially  for  steady-state  solutions,  they  must  be  discretized  both  temporally  and  spatially 
for  time  accurate  solutions.  First-order  discretization  gives  an  accurate  solution  for  cases 
with  simple  physical  and  flow  geometries,  but  for  most  cases  second-order  accurate 
solutions  are  more  desirable  (or  required).  Second-order  discretization  requires  values  at 
previous,  current  and  next  time  steps  for  an  unstructured  mesh.  The  under-relaxation 
parameter  scales  or  controls  updated  values  for  each  iteration,  thereby  stabilizing  the 
solution  but  slowing  convergence.  An  implicit  scheme  (where  all  flow  variables  are 
solved  simultaneously)  and  coupled  equations  lead  to  convergence  in  fewer  iterations 
than  an  explicit  scheme  and  segregated  equations,  respectively,  but  require  greater 
computational  resources  and  cannot  normally  be  run  in  a  parallel  computing  fashion. 
Current  state-of-the-art  flow  solvers  use  dual  time-stepping,  implicit-explicit  schemes  to 
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take  advantage  of  parallel  computing.  (FLUENT,  2001:  22.2.8,  22.4,  22.7.1,  22.9; 
Hoffman,  Chiang,  2000.11:  290) 

CFD  solutions  are  highly  mesh-dependent.  A  tradeoff  must  be  made  between 
overall  computation  time  and  solution  accuracy;  an  extremely  fine  mesh  will  typically 
give  an  accurate  solution  but  will  be  computationally  time-intensive,  while  a  coarse  mesh 
computes  quickly  but  may  give  an  inadequate  or  inaccurate  solution.  An  intermediate 
approach  is  to  use  a  coarse  mesh  which  may  be  adapted  by  redistributing  nodes  or  by 
adding  more  nodes  to  a  particular  region  for  better  resolution.  For  this  case,  FLUENT’S 
flow  software  provides  the  ability  to  use  a  roughly  converged  solution  to  determine 
where  best  to  refine  or  coarsen  the  mesh,  based  on  appropriate  flow  gradients,  distance 
from  the  wall,  user-specified  volumetric  regions,  etc.  The  flow  solver  then  continues 
toward  convergence  and  a  more  accurate  solution  with  this  modified  or  adapted  mesh, 
thus  optimizing  solution  accuracy  and  computational  efficiency  (FLUENT,  2001:  23.1). 

Turbulence  Models.  Turbulence  is  an  age-old  problem  that  is  relatively  easily 
observed,  as  noted  by  wise  King  Solomon,  but  not  easily  modeled.  He  wrote,  “The  wind 
goeth  toward  the  south,  and  turneth  about  unto  the  north;  it  whirleth  about  continually, 
and  the  wind  returneth  again  according  to  his  circuits”  ( The  Holy  Bible,  1979: 
Ecclesiastes  1:6). 

While  direct  numerical  simulation  (DNS)  best  predicts  turbulence  effects,  its  very 
fine  mesh  requirements  remain  too  computationally  demanding  -  given  today’s 
computing  capabilities  -  for  problems  with  complex  flow.  The  onus  is  on  the  user  to 
select  a  turbulence  model  which  gives  the  most  accurate  approximation  for  the  flow 
physics  in  a  specific  application.  Thus  for  this  application,  Ekaterinaris  and  Schiff 
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observed  that  “bursting  point  location  and  the  extent  of  the  vortex  breakdown  region  are 
seen  to  be  sensitive  to  the  turbulence  modeling”  (61).  As  part  of  this  study,  numerous 
turbulence  models  were  investigated.  These  models  include  Spalart-Allmaras  (S-A), 
Renormalization  Group  (RNG)  k-s,  Reynolds  Stress  (RSM),  and  Large-Eddy  Simulation 
(LES),  as  well  as  comparison  with  inviscid  and  laminar  flow  solutions.  It  must  be  noted 
that  in  general  these  turbulence  models  have  been  validated  for  cases  with  simple  flow 
features  and  simple  model  geometry  (Krai,  1998:  484);  while  this  case  involves  simple 
model  geometry,  the  flow  field  is  complex  so  careful  attention  must  be  afforded  in 
selecting  an  appropriate  model. 

While  it  is  included  for  comparison  and  has  shown  to  predict  a  relatively  accurate 
solution,  the  laminar  case  underestimates  the  strength  of  the  primary  and  secondary 
vortices  and  certainly  inaccurately  models  the  region  of  high  turbulence  within  and 
following  vortex  breakdown  (Ekaterinaris,  Schiff,  1990:  61;  Murayama,  Nakahashi, 
Sawada,  2001:  1311). 

Designed  for  aerospace  applications,  the  one-equation  S-A  turbulence  model 
solves  the  transport  equation  for  turbulent  or  kinematic  eddy  viscosity.  It  predicts 
solutions  best  for  cases  involving  wall-bounded  flow  and  adverse  pressure  gradients,  both 
of  which  are  part  of  this  study.  The  S-A  model  has  great  computational  efficiency  and 
has  been  shown  to  converge  on  a  solution  more  quickly  even  than  algebraic  turbulence 
models  which  encounter  discontinuities,  because  it  predicts  a  continuous  turbulence 
viscosity  distribution  (Mani,  Willhite,  Ladd,  1995;  Krai,  1998:  535).  While  it  is  also  well 
suited  for  meshes  with  unstructured  or  hybrid  boundary  layer,  the  S-A  model  was 
designed  for  problems  with  low  Reynolds  number,  attached  flow  and  mild  turbulence  and 
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is  still  relatively  new  in  its  validation  and  verification  (FLUENT,  2001:  10.2.4;  Spalart, 
Allmaras,  1992:1,  15).  However,  it  becomes  “very  complicated  and  often  ambiguous” 
when  attempting  to  model  physics  of  turbulent  flow  in  a  flow  field  with  significant  flow 
separation  and  unsteady  behavior  (Hoffman,  Chiang,  2000. Ill:  54). 

Other  numerical  investigations  have  shown  that  the  S-A  turbulence  model 
accurately  predicts  flow  properties  for  external  subsonic,  transonic  and  supersonic 
attached  flow,  and  for  supersonic  nozzle  and  impinging  jets  flow.  It  predicted  a  solution 
more  accurate  than  the  two-equation  k-s  model  for  supersonic  flow  over  a  flat  plate. 
While  the  S-A  model  has  successfully  predicted  qualitative  flow  physics,  its  quantitative 
predictions  are  considerably  less  accurate,  and  it  poorly  models  flow  in  turbulent  wake 
regions,  in  boundary  layers,  and  over  a  backward-facing  step  or  flat  trailing  edge. 
(Snyder,  Spall,  2000;  Krai,  1998;  Mani,  Willhite,  Ladd,  1995;  Spalart,  Allmaras,  1992: 
19) 

Accepted  nearly  industry-wide  as  an  accurate  turbulence  model  because  of  its 
strong  empirical  basis,  the  two-equation  k-s  model  (k  is  turbulent  kinetic  energy,  modeled 
theoretically;  e  is  turbulent  dissipation  rate,  modeled  empirically)  solves  two  independent 
transport  equations  and  accounts  for  compressibility  effects.  This  model  is  valid  only  for 
completely  turbulent  flow,  not  one  with  laminar  flow  in  various  regions,  but  the  RNG 
modification  reportedly  caters  to  time-dependent  turbulent  vortex  shedding  (FLUENT, 
2001:  10.2.5,  10.2.11,  10.4.1).  Numerical  studies  have  shown  that  the  k-s  model  more 
accurately  predicts  flow  in  a  compressible,  turbulent  field  than  does  the  S-A  model  and 
that  it  is  well  suited  to  cases  with  supersonic  attached  flow,  internal  flow  and  external 
turbulent  boundary  layers.  Disadvantages  include  poor  accuracy  for  cases  with  vortical 
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and  rotating  flow  and  with  flow  in  non-circular  ducts  and  turbulent  wakes;  also  the  k-s 
model  requires  about  10%  more  computational  memory  and  2-7  times  more  processing 
time  than  does  the  S-A  model.  (Krai,  1998;  Krai,  Mani,  Ladd,  1996;  Versteeg, 
Malalasekera,  1995:  75) 

Similar  to  k-s  is  the  k-co  turbulence  model,  which  solves  the  transport  equations 
for  k  and  to,  or  specific  dissipation  rate.  FLUENT  provides  two  versions  of  the  k-co 
model:  the  Standard  k-co  model  accounts  for  effects  from  compressibility  and  low 
Reynolds  number;  the  Shear-Stress  Transport  (SST)  k-co  model  accounts  for  the  principal 
turbulent  shear  stress  and  combines  benefits  of  both  k-co  and  k-s  models  within  the 
boundary  layer  (FLUENT,  2001:  10.5,  10.5.1-2).  In  computational  tenns,  both  k-co 
models  require  memory  comparable  to  that  required  by  the  k-s  model,  but  the  SST  k-co 
model’s  processing  time  per  iteration  is  on  the  order  of  that  required  for  S-A 
computations  -  due  to  its  limiter  on  turbulent  viscosity.  Both  k-co  models  predict 
relatively  accurately  supersonic  external  flow  with  an  attached  boundary  layer;  the  SST 
version  does  well  also  with  supersonic  separated  flow  and  nozzle  flow,  but  in  these  cases 
it  does  not  predict  better  than  does  the  S-A  model.  The  SST  k-co  model  poorly  predicts 
flow  that  is  subsonic  and  highly  separated  (Krai,  1998:  484,  535-538).  Since  the  k-co 
models  are  only  a  small  improvement  over  the  k-s  model,  largely  comparable  to  the  S-A 
model,  and  inferior  to  the  RSM  model,  they  were  not  evaluated  in  this  research  effort. 

More  complex  than  the  other  models,  the  seven-equation  RSM  solves  seven 
transport  equations  -  five  for  the  Reynolds  stresses  and  two  for  k  and  s.  It  accounts  for 
flow  rotation,  compressibility,  curved  streamlines,  the  anisotropic  or  directional 
characteristic  of  turbulence  in  swirling  flows,  and  it  does  not  assume  that  “turbulent 


1-22 


stresses  respond  immediately  to  changes  in  the  mean  strain  rate”  (Krai,  1998:  483). 
However,  it  is  limited  by  potential  inaccuracies  of  the  empirical  e  model,  assumes  local 
homogeneity  and  equilibrium  in  turbulence,  and  requires  50-60%  more  computation  time 
per  iteration  and  15-20%  more  memory  than  the  other  turbulence  models  (FLUENT, 
2001:  10.2.10-11;  Krai,  1998;  Versteeg,  Malalasekera,  1995:  78).  For  steady, 
incompressible  flow  at  Re  =  8.0  x  1 05,  the  RSM  model  was  verified  as  more 
quantitatively  accurate  than  the  S-A  model  for  flow  within  the  boundary  layer  and  the 
turbulent  vortex  structure  (Snyder,  Spall,  2000). 

While  the  S-A,  k-s  and  RSM  turbulence  models  use  Reynolds-Averaged  Navier- 
Stokes  equations,  which  model  mean  flow  quantities  and  which  generally  require  less 
computational  time  (for  steady  and  unsteady  flows)  than  DNS,  LES  uses  an  alternative 
filtering  approach  to  solving  for  the  flow  variables  (unsteady  flow  only).  By 
appropriately  modifying  the  Navier-Stokes  equations,  turbulent  eddies  smaller  than  the 
filter  (or  grid  cell  size)  are  removed  and  approximated  with  an  isotropic  model,  while  the 
large  eddies  are  directly  resolved  through  the  discretized  Navier-Stokes  equations.  Small 
eddies  are  more  universal  and  depend  less  on  geometery  and  are  thus  more  easily 
modeled,  whereas  large  eddies  are  specific  to  a  problem’s  geometry  and  boundary  and 
initial  conditions.  Disadvantages  of  LES  include:  it  requires  a  fine  numerical  mesh, 
which  results  in  high  computational  cost;  along  the  walls,  the  mesh  must  be  particularly 
fine,  since  LES  basically  uses  direct  numerical  simulation  in  this  region;  it  has  not  been 
well  tested  for  cases  with  other  than  simple  geometries;  and  it  does  not  account  for 
compressible  flow,  though  it  allows  for  variable  density.  This  study  used  the 
Smagorinsky-Lilly  Model,  which  is  accurate  in  many  cases  for  flow  with  high  Reynolds 
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number  (that  is,  flow  in  regions  other  than  near  walls)  to  model  the  small  or  subgrid 
turbulent  eddies.  (FLUENT,  2001:  10.2.1,  10.7,  10.7.2;  Mathieu,  Scott,  2000:  340-353) 
Solution  Analysis  and  Evaluation 

Once  convergence  or  fully  developed  flow  is  realized,  solution  data  may  be 
analyzed  and  evaluated  using  any  number  of  visualization  software  packages.  In  this 
case,  variables  of  interest  include  pressure,  velocity,  vorticity  (rotation  measurement), 
helicity  (dot  product  of  vorticity  and  velocity  vector),  and  turbulence.  Flow  visualization 
with  such  software  may  be  accomplished  virtually  instantaneously,  whereas  wind  tunnel 
flow  visualization  requires  injection  of  smoke  or  dies  and  use  of  lasers,  high-speed 
cameras,  etc.,  and  consumes  a  good  deal  of  time  and  resources. 

If  convergence  has  not  been  achieved  or  if  the  solution  is  otherwise  deemed 
inaccurate,  possible  actions  include  executing  more  iterations  with  smaller  required 
residuals,  decreasing  the  CFL  number  or  time  step,  decreasing  the  under-relaxation 
parameter(s),  feeding  a  converged  first-order  solution  into  a  second-order  computation, 
refining  the  mesh,  and/or  selecting  a  different  turbulence  model. 

Previous  Work  on  Numerical  Simulation  of  Vortex  Breakdown  Control 
Other  numerical  studies  of  vortex  breakdown  typically  assume  axisymmetric, 
laminar,  incompressible,  steady  flow  at  low  Reynolds  number  (Murayama,  Nakahashi, 
Sawada,  2001;  Novak,  Sarpkaya,  1999:  825;  Ekaterinaris,  Schiff,  1990;  Leibovich,  1978: 
243;  and  numerous  studies  described  and  cited  in  these  sources).  While  some  studies 
have  investigated  unsteady  flow,  one  expert  who  has  researched  and  investigated  vortex 
breakdown  for  more  than  30  years,  stated,  “...there  has  not  yet  been  a  turbulence  model 
capable  of  dealing  with  nonisotropic  turbulence  in  swirling  flows  ...  subjected  to 
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streamline  curvature  and  strong  radial  pressure  gradients”  (Novak,  Sarpkaya,  1999:  825, 
834). 

While  these  studies  have,  in  general,  closely  matched  numerical  and  experimental 
results,  no  numerical  study  has  investigated  the  specific  set  of  conditions,  particularly 
with  along-core  blowing,  applied  to  this  specific  geometry.  Also,  the  studies  referenced 
above  modeled  the  delta  wing  in  farfield,  freestream  conditions,  as  opposed  to  including 
the  wind  tunnel  geometry  as  part  of  the  model. 
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II.  Delta  Wing  Model  and  Computational  Facilities 


The  scope  of  this  branch  of  the  overall  research  effort  does  not  include  wind 
tunnel  test  article  data  generation  and  collection,  but  the  physical  configuration  at  UC  is 
described  hereafter  to  show  correlation  with  the  numerical  representation. 

Delta  Wing  Model  and  Facilities 

The  model,  pictured  in  Figure  2.1,  is  a  half-span  (port  side),  aluminum  delta  wing 
with  sharp  leading  edge,  sweep  angle  of  60  degrees,  flat  trailing  edge,  and  no  camber.  It 
is  1.27  cm  thick  with  30-degree  leading  edge  bevel  and  has  a  removable  top  plate  and 
hollow  interior  to  accommodate  pressure  sensor  instrumentation  and  tubing  for  the 
blowing  ports.  Root  chord  length  (c)  measures  34.3  cm  and  half-span  (s,  where  span,  b  = 


Figure  2.1-  Image  of  Half  Delta  Wing  Test  Article 
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2s)  is  19.8  cm.  Wing  upper  surface  (or  planform)  area  is  339.57  cm'  and  aspect  ratio  is 
2.31.  The  leading  edge  is  not  perfectly  sharp,  as  it  appears  to  have  been  ground  slightly. 
Further,  the  wing  is  many  years  old  and  thus  has  numerous  small  abrasions,  nicks  and 
dents  along  the  leading  edge  and  over  the  entire  surface. 

The  three  blowing  ports,  located  chordwise  at  x/c  =  0.30,  0.60  and  0.80,  and 
spanwise  at  y/s  =  0.21,  0.41  and  0.54  were  positioned  based  on  visual  placement  under 
the  observed  vortex  core  centerline.  Each  blowing  port  has  a  straight  nozzle  drilled  in  a 
circular  disk,  where  nozzle  dimensions  are  0.9525  cm  (3/8  in.)  in  length  and  0.0794  cm 
(1/32  in.)  in  diameter.  Based  on  results  from  the  study  at  Louisiana  State  University, 
using  the  same  wing,  relative  optimum  blowing  angles  were  established  for  pitch  and 
azimuthal  directions,  where  pitch  angle  is  35  degrees,  measured  from  the  wing  surface, 
and  azimuthal  angle  is  155  degrees,  measured  counter-clockwise  (positive)  from  a  line 
parallel  to  the  wing’s  root,  as  shown  in  Figure  2.2  (Guillot,  1999:  21-22).  Each  blowing 
port  is  connected,  via  6  m  of  rubber  hose,  to  a  source  of  compressed  air  at  available 
values  of  586,  483,  345  or  207  kPa  (85,  70,  50  or  30  psig,  respectively). 


Figure  2.2  -  Jet  Blowing  Angle:  a)  Azimuthal,  b)  Pitch  Directions  (Guillot,  1999:  21) 
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Thirty  pressure  taps  on  the  wing  surface  are  located  in  straight  spanwise  lines  and 
in  chordwise  positions  along  x/c  =  0.35,  0.55,  0.75  and  0.95.  Tubes  from  each  pressure 
tap  and  blowing  port  exit  the  test  article  through  a  mounting  bracket  3.81  cm  (1.5  in.) 
long  and  1.27  cm  (0.5  in.)  thick,  with  rounded  front  and  back  ends  of  radius  0.635  cm 
(0.25  in.).  Static  pressure  measurement  error  on  the  wing  surface  is  on  the  order  of  a  few 
Pascal,  based  on  0.8  to  0.9  Pa  instrument  error  (maybe  more  since  the  transducers  are 
about  two  decades  old),  in  addition  to  tube  losses.  Each  pressure  measurement  was  time 
averaged  from  5  seconds  of  data  collected  at  100  Hz  sampling  rate  (May,  2002. a). 

The  wing  is  mounted  at  an  angle  of  attack  equal  to  approximately  15  degrees  and 
flush  to  a  boundary  layer  refreshing  plate,  0.635  cm  thick,  45.72  cm  (18  in.)  long,  and 
30.48  cm  (12  in.)  high,  shown  in  Figures  2.3  and  2.4.  This  plate  has  a  45-degree  bevel 
and  sharp  leading  edge,  where  the  wing  is  mounted  to  the  underside;  the  wing  apex  is 
located  about  5  cm  downstream  of  the  refresher  plate’s  leading  edge.  The  plate  is  offset 
from  the  wind-tunnel  wall  by  1 .27  cm,  which  effectively  compresses  into  this  gap  the 


Figure  2.3  -  Half  Delta  Wing  with  Boundary  Layer  Refresher  Plate 
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Figure  2.4  -  Partial  Isometric  View  of  Wing  and  Boundary  Layer  Refresher  Plate 
wind-tunnel  wall  boundary  layer  and  gives  laminar  freestream  conditions  at  the  wing 
apex.  This  serves  as  an  approximation  of  freestream  conditions  in  unconfined  flow. 

Wind  tunnel  test  section  dimensions  are  60.96  cm  (24  in.)  high  by  60.96  cm  wide. 
The  test  article  is  mounted  about  200  cm  downstream  of  the  relatively  steady  and  evenly 
distributed  flow  inlet  and  about  300  cm  upstream  of  the  wind  tunnel  exit/recirculation 
chamber.  Figure  2.5  shows  tunnel  and  boundary  layer  refresher  plate  dimensions. 


Figure  2.5  -  Boundary  Layer  Refresher  Plate  and  Wind  Tunnel  Test  Section 
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Test  data  were  generated  in  a  UC  lab,  where  ambient  pressure  was  measured  at 
99.56  kPa  and  ambient  temperature  was  298  K.  Dynamic  pressure  measurements  inside 
the  wind  tunnel  and  upstream  of  the  wing  averaged  about  143  Pa  (May,  2002. a).  These 
pressure  measurements  correspond  to  an  air  density  of  1.208  kg/m  and  freestream 
velocity  of  15.4  m/s  (U.S.  Government  Printing  Office,  1976),  which  then  gives  a 
Reynolds  number  of  3.56  x  1 05,  based  on  root  chord  length,  and  a  total  temperature  of 
298. 12  K  inside  the  wind  tunnel.  Refer  to  Appendix  A  for  raw  data. 

Computational  Facilities 

Computational  research  and  investigation  for  this  study  was  performed  in  the 
Computational  Fluid  Dynamics  Laboratory  (CFD  Lab)  of  the  Department  of  Aeronautics 
and  Astronautics  and  in  the  UNIX  Computer  Lab,  both  of  the  Graduate  School  of 
Engineering  and  Management,  Air  Force  Institute  of  Technology,  using  the  following 
hardware  resources  and  commercially  available  software. 

Hardware 

Mesh  generation,  early  flow  solutions  and  most  post-processing  were  performed 
on  Dell  Precision  530  workstations,  using  Redhat  Linux  7.2  or  8.0  operating  systems. 
Two  of  the  three  available  workstations  in  the  CFD  Lab  feature  dual  1 .4  GHz  Pentium  4 
Xeon  processors,  512  MB  of  RAM  and  37  GB  of  hard  disk  space,  and  the  third  features 
dual  2.3  GHz  Pentium  processors,  2,048  MB  of  RAM  and  over  120  GB  of  disk  space. 
They  are  connected  to  a  16-node  Beowulf  processing  cluster  (also  known  as  the  Hydra 
cluster)  by  Aspen  Systems,  Inc.  (http://www.aspsys.com),  where  each  node  has  two  1.2 
GHz  Athlon  processors  and  1,024  MB  of  RAM.  Computing  was  also  performed  on  a  64- 


2-5 


node  Beowulf  cluster  (known  as  the  Aspen  cluster),  where  each  node  has  two  1.0  GHz 
processors  and  1,024  MB  of  RAM.  The  Aspen  cluster  was  accessed  remotely  from  a  Sun 
Ultra  80  UNIX  workstation,  using  a  Sun  operating  system  and  featuring  2,048  MB  of 
RAM  and  18  GB  of  hard  drive  space.  Additional  post-processing  was  performed  on  a 
Dell  Precision  530  workstation,  using  a  Microsoft  Windows  2000  Professional  operating 
system  and  featuring  a  1.5  GHz  Pentium  4  processor  and  512  MB  of  RAM. 

Software 

Commercially  available  software  was  used  for  mesh  generation,  flow 
initialization  and  computation,  and  post-processing.  The  hybrid  mesh,  or  combination  of 
structured  and  unstructured  grids,  was  generated  using  Gridgen  Version  14.02 
(Pointwise,  Inc.,  http://www.pointwise.com).  All  flow  computations  were  performed 
with  FLUENT  Version  6.0.20,  which  is  a  general-purpose  flow  solver  (FLUENT,  Inc., 
http://www.fluent.com).  Since  it  uses  a  finite-volume  discretization  scheme,  FLUENT  is 
well  suited  for  solving  unstructured  or  hybrid  domains.  FLUENT  has  numerous  options 
for  computational  schemes,  including  steady  or  unsteady  flow,  first-  and  second-order 
spatial  and  temporal  accuracy,  incompressible  or  compressible  flow,  coupled  or 
segregated  equation  solver,  explicit  or  implicit  numerical  scheme,  and  inviscid  or  viscous 
flow,  featuring  each  of  the  turbulence  models  discussed  in  Chapter  1  (FLUENT,  2001). 
Results  from  and/or  justification  for  selected  options  in  FLUENT  follow  in  Chapter  3. 
Post-processing  was  done  both  with  FLUENT  and  with  FIELD  VIEW  Version  8.0 
(Intelligent  Light,  http://www.ilight.com).  FLUENT  has  an  option  to  export  solution  data 
in  FIELD  VIEW  format.  Post-processing  included  generating  plots,  visualizing  isometric 
contours  of  various  flow  variables,  generating  multiple  two-dimensional  slices  in  one 
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image,  creating  streamlines,  surface  flow  and  surface  pressure  contours,  and  identifying 
vortex  core  centerline. 
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III.  Numerical  Simulation  without  Flow  Control 


As  discussed  in  Chapter  I,  tackling  a  CFD  problem  consists  of  three  iterative  steps 
-  mesh  generation,  numerical  modeling  and  computation,  and  solution  analysis  and 
evaluation.  As  the  process  is  iterative,  it  also  involves  mesh  refinement  and/or  adaptation 
and  possible  variation  of  solver  parameters,  flow  models,  and  initial  and  boundary 
conditions.  Discussion  and  comparisons  follow,  to  include  parallel  computing.  Chapters 
III  and  IV  address  these  steps  for  this  problem,  where  Chapter  III  discusses  the  baseline 
case  (half  delta  wing  with  no  along-core  blowing)  and  Chapter  IV  presents  an  along-core 
blowing  case. 

Mesh  Generation 

Since  the  half  delta  wing  test  article  is  a  number  of  years  old  and  was  passed  from 
one  university  to  another,  the  original  technical  drawing  was  no  longer  available  resulting 
in  possibly  inaccurate  measurements  of  wing  dimensions  and  blowing  port  locations. 
Furthermore,  the  wind  tunnel  setup  was  no  longer  fully  assembled  or  being  used  during 
the  time  of  initial  mesh  development,  such  that  some  details  of  the  setup  were  initially 
absent,  including  existence  of  the  boundary  layer  refreshing  plate  and  precise  angle  of 
attack  and  freestream  velocity  measurements.  Thus  the  initial  mesh  and  its  first  10 
revisions  all  lack  inclusion  of  a  refreshed  boundary  layer  about  5  cm  upstream  of  the 
wing  apex.  However,  some  experience  was  gained  in  passing  through  the  various  mesh 
iterations,  so  they  are  included  here  for  pedagogical  edification.  Mesh  details,  including 
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number  of  node  points,  faces  and  cells,  are  discussed  qualitatively  (or  approximately) 
until  discussion  of  the  final  configuration. 

While  structured  grids  traditionally  give  more  accurate  solutions  for  relatively 
simple  geometries,  unstructured  meshes  may  achieve  a  sufficiently  accurate  solution  with 
fewer  cells,  less  construction  effort  and  hence  less  computation  time.  After  significant 
effort  to  create  a  structured  grid,  to  include  meshing  discontinuities  caused  by 
introduction  of  blowing  ports,  it  was  determined  an  unstructured  approach  would  be 
quicker,  more  easily  revised  as  needed  and  would  likely  provide  acceptable  results. 

The  initial  mesh  consisted  of  the  half  delta  wing  placed  flush  against  a  wall  which 
extended  one  chord  length  downstream  of  the  trailing  edge,  one  chord  length  upstream  of 
the  apex,  and  one  chord  length  both  above  and  below  the  wing  in  a  horseshoe-shaped,  C- 
type  grid.  It  was  assumed  unnecessary  to  model  the  wind  tunnel  walls,  rather  model  the 
wing  in  freestream  conditions.  For  this  and  with  ensuing  discussion,  a  chord  length  is 
defined  as  25  cm,  whereas  the  root  chord  is  34.3  cm.  The  domain  extended  spanwise 
about  7  cm  beyond  the  wingtip,  where  leading  and  trailing  edges  meet.  All  nodes  were 
equally  and  relatively  densely  spaced,  such  that  the  volume  consisted  of  about  300,000 
tetrahedral  cells  and  60,000  nodes.  The  wing,  which  remained  virtually  unchanged 
throughout  later  grid  modifications,  had  the  domain’s  largest  concentration  of  triangular 
faces,  to  better  capture  boundary  layer  effects  and  to  have  higher  resolution  on  the 
surface  where  experimental  pressure  sensors  were  located.  Angle  of  attack  was 
simulated  by  setting  Cartesian  components  for  the  freestream  velocity.  Thus  a  new  mesh 
would  not  need  to  be  generated  for  variations  in  angle  of  attack,  rather  the  initial 
conditions  would  be  appropriately  modified. 
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The  origin  was  located  on  the  root  chord,  incident  with  the  upper  wing  surface, 
and  at  x/c  =  0.55,  which  is  aligned  with  one  of  the  rows  of  pressure  sensors  on  the  test 
article.  Axes  were  oriented  such  that  the  positive  x-axis  followed  the  direction  of  air 
flow,  positive  v-axis  traveled  outboard  of  the  origin,  and  positive  z-axis  pointed 
downward  from  or  nonnal  to  the  wing’s  lower  surface.  This  convention  was  maintained 
throughout  all  later  grid  revisions. 

Impetus  for  modifying  the  grid  came  from  the  following:  comparing  plots  of  the 
numerically  predicted  pressure  coefficients  (Cp)  with  those  detennined  experimentally  at 
the  30  locations  described  in  Chapter  II;  viewing  surface  (two-dimensional,  as  well  as 
isometric)  contours  of  flow  variables  along  the  wall  and  within  the  domain  to  ensure  the 
solution  contours  were  captured  and  resolved  within  the  domain  extents;  determining 
whether  the  solution  was  sufficiently  converged  or  flow  was  fully  developed  (via  error 
residuals  for  continuity,  velocity,  energy  and  viscosity  equations,  via  mass  flow 
conservation,  and  via  constant  (for  steady  cases)  or  steady  cyclic  (for  unsteady  cases)  CL 
and  CD);  and  comparing  CL  with  experimental  and  theoretical  values.  Such  plots, 
contours  and  other  post-processing  graphics  are  compared  and  evaluated  later  in  this 
chapter  to  give  justification  for  grid  revisions;  thus  detailed  reasoning  for  grid 
modification  is  not  presented  in  this  section. 

Another  judge  of  numerical  mesh  suitability  is  to  compute  y+  (or  y-plus)  values 
along  the  walls  when  using  turbulence  models.  Y-plus  indicates  whether  there  is 
adequate  grid  resolution  near  the  wall  and  is  defined  by 
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where  aT  =  (xjpw)  ~  is  friction  velocity,  zw  is  wall  shear  stress,  pw  is  fluid  density  at  the 
wall,  yp  is  distance  from  point  P  to  wall,  p  is  fluid  density  at  point  P,  and  p  is  fluid 
viscosity  at  point  P.  The  Spalart-Allmaras  (S-A),  k-s  and  Reynold  Stress  (RSM) 
turbulence  models  have  a  desirable  wall  y+  range  of  30  <y+  <  60,  while  the  Large  Eddy 
Simulation  (LES)  model  requires  y+  ~  1  along  the  walls  (FLUENT,  2001:  10.8.1,  10.9.1- 
5,  27.4;  Hoffman,  Chiang,  2000. Ill:  55). 

Revisions  are  denoted  alphabetically,  so  Revision  A,  shown  in  Figure  3.1, 
extended  the  domain  20  chord  lengths  downstream  and  five  chord  lengths  above  and 


Figure  3.1  -  Mesh  Generation,  Revision  A:  Farfield  Model  with  One  Tunnel  Wall 
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below  the  wing,  but  the  upper  and  lower  trailing  extents  were  connected  via  parabola  vice 
the  horseshoe  shape  of  the  initial  mesh.  Node  distribution  on  the  wing  remained 
unchanged,  though  the  wing  was  offset  1.27  cm  from  the  wall  to  account  for  the 
mounting  bracket  and  to  move  the  wing  farther  from  the  wall  boundary  layer.  The  rest  of 
the  domain  had  considerably  fewer  cells  and  nodes,  reducing  to  about  80,000  tetrahedral 
cells  and  20,000  nodes.  In  this  revision,  nodes  were  clustered  toward  the  leading  end  (or 
parabola  base)  to  give  better  resolution. 

Revision  B,  shown  in  Figure  3.2,  extended  downstream  to  100  chord  lengths  past 
the  wing  trailing  edge  and  up  an  additional  20  chord  lengths  beyond  the  wing  upper 
surface  -  basically  enlarging  the  computational  domain  to  capture  the  generated  vortex. 
About  58,000  nodes  were  added,  mostly  in  the  region  around  the  wing  as  seen  in  Figure 
3.3,  to  further  enhance  resolution,  resulting  in  a  total  of  about  400,000  cells. 


Figure  3.2  -  Mesh  Generation,  Revision  B:  Upper  and  Downstream  Extensions 
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Figure  3.3  -  Mesh  Generation,  Revision  B:  Closer  View  of  Region  around  Wing 
Revision  C  expanded  the  downstream  extent  to  150  chord  lengths  (37.5  m)  and 


upper  extent  to  100  chord  lengths  (25  m),  as  seen  in  Figure  3.4.  Number  of  nodes  was 
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basically  maintained  in  the  region  around  the  wing  but  greatly  reduced  in  the  farfield 
areas,  bringing  total  cell  count  to  246,000,  with  48,000  nodes.  Note  the  relatively 
minuscule  size  of  the  wing  in  this  domain. 

Revision  D,  shown  in  Figure  3.5,  simply  extended  the  spanwise  dimension  of  the 
domain  to  60  chord  lengths  and  nearly  maintained  the  same  number  of  nodes  and  cells. 
The  parabolic  front  end  was  only  continued  for  a  few  chord  lengths  in  the  spanwise 
direction  before  it  flattened,  which  is  characteristic  of  the  automated  grid  generation;  to 
obviate  this  problem,  the  user  would  need  to  generate  a  more  closely  ribbed  structure 
(which  was  done  in  Revision  F)  or  integrate  a  structured  grid  into  this  portion  of  the 
domain.  This  domain  size  effectively  “captured”  all  solution  contours,  but  computed  data 
still  did  not  adequately  match  experimental  results. 


Figure  3.5  -  Mesh  Generation,  Revision  D:  Spanwise  Extension 
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Meshes  to  this  point  had  not  really  accounted  for  flow  disturbances  which  may 
propagate  upstream  of  the  wing,  since  flow  was  low  subsonic  in  this  case.  Thus  Revision 
E  extended  the  domain  20  chord  lengths  upstream  of  the  wing,  added  another  10  chord 
lengths  to  the  extent  below  the  wing,  as  shown  in  Figure  3.6,  and  dramatically  increased 
the  number  of  nodes  in  the  wing  region,  as  shown  in  Figure  3.7.  Further,  the  wing  was 
separated  as  an  entity  to  allow  for  calculation  of  Cl  and  Cd .  These  changes  brought  the 
number  of  tetrahedral  cells  to  523,000  and  number  of  nodes  to  100,000. 


Figure  3.6  -  Mesh  Generation,  Revision  E:  Upstream  and  Fower  Extensions 
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Figure  3.7  -  Mesh  Generation,  Revision  E:  Closer  View  of  Region  around  Wing 


Since  Revision  E  gave  an  unacceptable  solution,  Revision  F  extended  the  domain 
downstream  to  250  chord  lengths  but  reduced  the  number  of  tetrahedral  cells  to  448,000, 
with  88,000  nodes.  This  revision  also  eliminated  possible  inaccuracy  from  the  flattened 
front  side  by  making  it  parabolic  to  the  spanwise  extent,  as  seen  in  Figure  3.8.  These 
changes  had  little  effect  on  the  solution,  since  the  domain  was  already  sufficiently  large 
and  because  flow  far  outboard  of  the  wing  had  no  recirculation  issues. 
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Figure  3.8  -  Mesh  Generation,  Revision  F:  Parabolic  Upstream  Extent 


At  this  point,  it  was  determined  the  most  prudent  course  of  action  would  be  to 
return  to  the  drawing  board  and  create  the  domain  differently.  That  is,  modeling  the 
computational  domain  to  be  the  same  as  the  experiment’s  physical  domain  -  a  wind 
tunnel  -  was  the  basis  for  Revision  G. 

Revision  G  eliminated  all  existing  boundaries  around  the  wing  and  created  a  four¬ 
sided  wind  tunnel  test  section  60.96  cm  in  height,  60.96  cm  in  width,  and  which  extended 
20  chord  lengths  both  upstream  and  downstream  of  the  wing.  Still  offset  1 .27  cm  from 
the  tunnel  wall,  the  wing  was  rotated  15°  about  the  y-axis  to  simulate  angle  of  attack, 
since  air  would  need  to  flow  straight  through  the  tunnel  or  perpendicular  to  the  inlet  and 
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outlet  boundaries.  A  problem  with  rotating  the  wing  to  the  proper  angle  of  attack  is  that 
any  modification  of  that  angle  requires  a  new  grid  to  be  generated,  as  opposed  to  simply 
altering  a  boundary  condition  in  the  flow  solver.  This  setup  called  for  a  hybrid  mesh, 
consisting  of  two  structured  blocks  and  one  unstructured  block.  The  unstructured  block 
centered  around  the  wing,  extending  two  chord  lengths  upstream  and  downstream  of  the 
origin  and  filling  that  section  of  tunnel,  as  shown  in  Figure  3.9.  Structured  blocks 
extended  outward  from  each  end  of  the  unstructured  volume,  and  nodes  were  distributed 
more  densely  within  a  region  5  cm  from  each  wall  surface  to  better  capture  boundary 
layer  effects.  The  new  domain  consisted  of  about  100,000  hexahedral  cells  (structured 


Figure  3.9  -  Mesh  Generation,  Revision  G:  Hybrid  Grid  Including  Wind  Tunnel  Model 
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volume)  and  300,000  tetrahedral  and  pyramid  cells  within  and  along  the  borders  of  the 
unstructured  volume. 

Revision  H  enhanced  wall  boundary  layer  resolution  within  the  unstructured 
block  by  isolating  the  unstructured  from  structured  blocks  via  “transition”  blocks  and  by 
creating  an  interior  volume  around  the  wing.  Referring  to  Figure  3.10,  note  the 
structured  walls  have  fewer  nodes  than  were  needed  to  adequately  resolve  boundary  layer 
effects  in  the  unstructured  wing  section  of  the  tunnel.  Therefore,  a  transition  zone  (about 
5  cm  long)  was  put  in  place,  where  the  node  count  was,  for  example,  3 1  on  one  side  and 


Figure  3.10  -  Mesh  Generation,  Revision  H:  Enhanced  Boundary  Layer  Resolution 
(Features  Some  Nodal  Dimensions  for  Revision  L) 


3-12 


7 1  on  the  other;  the  transition  zone  had  adequate  resolution  but  was  short  in  length  while 
the  unstructured  volume  then  had  good  resolution.  Next,  an  interior  volume  was  created 
5  cm  inside  all  unstructured  tunnel  walls,  shown  as  purple  mesh  in  Figure  3.10,  except  on 
the  wing-mounted  wall  where  resolution  was  already  adequate.  This  resulted  in  a 
boundary  layer  3-7  cells  thick  throughout  the  model’s  unstructured  portion.  Total  cell 
count  increased  to  590,000. 

In  an  effort  to  reduce  boundary  layer  turbulence  which  engulfed  the  wing  in 
Revision  H,  Revision  I  truncated  the  leading  structured  tunnel  to  one  chord  length,  such 
that  the  pressure  inlet  was  then  about  60  cm  from  the  wing  apex,  as  seen  in  Figure  3.11. 


At  this  point,  discussion  with  the  UC  test  engineer  revealed  the  presence  of  a 
boundary  layer  refresher  plate  during  wind  tunnel  testing  (May,  2002. c),  which  generated 
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essentially  laminar  flow  near  the  wing  apex,  as  discussed  in  Chapter  II.  For  Revision  J 
the  refresher  plate’s  effect  was  approximated  by  truncating  the  domain  to  7.6  cm 
upstream  -  later  corrected  by  UC  test  engineer  to  be  5  cm  (May,  2002.b)  -  of  the  wing 
apex  and  providing  evenly  distributed,  non-turbulent  flow  at  the  inlet,  as  seen  in  Figure 
3.12.  This  revision  also  added  more  nodes  to  the  unstructured  volume,  bringing  total 
number  of  cells  to  nearly  800,000.  The  wing  was  still  offset  from  the  wall,  an  error 
corrected  in  the  next  revision. 


Figure  3.12  -  Mesh  Generation,  Revision  J:  More  Truncated  Distance  to  Inlet 


Revision  K  placed  the  wing  flush  against  the  tunnel  wall/refresher  plate,  modified 


the  inlet  to  be  5  cm  upstream  of  the  wing  apex,  and  included  several  derivative  meshes  at 
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various  angles  of  attack.  Figure  3.13  shows  this  numerical  mesh  with  wing  at  15°  angle 
of  attack.  With  wing  against  wall,  the  number  of  domain  cells  reduced  to  764,000. 


Figure  3.13-  Mesh  Generation,  Revision  K:  Wing  Flush  Against  Wall 


From  Revision  K  findings,  it  was  determined  that  the  angle  of  attack  should  be 
1 8°  for  best  results;  more  detail  is  provided  later  in  this  chapter.  Revision  L  incorporated 
this  modified  wing  angle  and  increased  fidelity  of  the  boundary  layer  refresher  plate 
approximation  by  including  the  plate’s  leading  edge  45°  bevel,  coupled  with  the  re¬ 
addition  of  structured  upstream  tunnel,  as  seen  in  Figure  3.14.  About  4.3  cm  upstream  of 
the  wing  apex,  the  tunnel  wall  angles  45°  inboard  for  0.635  cm  in  the  x-direction,  then 
transitions  for  a  centimeter  or  so  into  a  structured  mesh,  which  then  continues  to  a 
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Figure  3.14  -  Mesh  Generation,  Revision  L:  Tunnel  with  Nodal  Dimensions 
distance  200  cm  upstream  of  the  origin.  The  tunnel  runs  500  cm  downstream  of  the 
origin,  which  is  200  cm  greater  than  the  actual  wind  tunnel  test  section’s  downstream 
segment,  but  this  was  done  to  allow  greater  extent  for  flow  resolution  and  accuracy. 

The  numerical  boundary  layer  refresher  plate  is  an  approximation  for  three 
reasons:  first,  this  ramp  extends  to  the  tunnel  ceiling  and  floor  vice  covering  only  half  the 
vertical  distance  of  the  tunnel  wall,  as  shown  in  Figure  2.5;  second,  this  simulated  plate  is 
not  offset  by  1.27  cm  from  the  tunnel  wall;  and  third,  the  approximated  plate  runs  the 
remaining  downstream  distance  of  the  tunnel,  instead  of  extending  only  45.72  cm  from 
the  plate’s  leading  edge.  The  first  and  third  approximations  likely  produce  negligible 
differences  in  the  solution,  but  approximating  the  plate  as  flush  against  the  tunnel  wall  is 
not  the  same  as  the  actual  wind  tunnel  setup  in  the  vicinity  of  the  wing.  Flere  a  tradeoff 
was  made  for  greater  modeling  simplicity,  as  the  solver  would  possibly  have  needed  to 
deal  with  additional  flow  physics  between  the  refresher  plate  and  tunnel  wall  and  with 
more  complicated  downstream  flow.  Revision  O  removed  the  second  and  third 
approximations. 
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Figure  3.14  shows  nodal  dimensions  for  the  tunnel  walls  and  inlet  and  outlet 
faces,  and  Figure  3.15  shows  them  for  the  transition  zone  between  structured  and 
unstructured  segments,  for  the  unstructured  volume  walls  and  for  the  wing.  To  avoid  a 


Figure  3.15  -  Mesh  Generation,  Revision  L:  Nodal  Dimensions  for  Tunnel  Wing  Section 
superfluous  figure,  nodal  dimensions  for  Revision  L’s  interior  volume  surrounding  the 
wing  are  shown  in  Figure  3.10,  since  this  grid  structure  remained  the  same;  faces  not 
shown  have  essentially  the  same  nodal  dimensions  as  the  bottom  face.  This  interior 
volume  allowed  for  a  5 -cm  wall  boundary  layer  thickness  (where  cells  were  more 
concentrated)  throughout  the  unstructured  zone  of  the  tunnel.  Furthermore,  in  the 
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structured  portions  of  the  tunnel,  the  model  provided  a  5 -cm  zone  along  the  tunnel  walls 
for  boundary  layer  resolution.  Length  of  that  zone  was  divided  into  7  square  cells,  where 
the  one  closest  to  the  wall  measured  0.2  cm  on  each  side,  then  they  increased  to  1.3  cm 
for  the  seventh  cell. 

Not  including  interior  faces,  the  numerical  model  consists  of  67,598  triangular 
(unstructured)  and  quadrilateral  (structured)  faces;  this  number  does  include  faces  of  the 
interior  volume  from  Figure  3.10.  The  wing  claims  24.2%  of  that  number  of  faces. 
There  are  a  total  of  585,297  cells  and  176,626  nodes;  the  cells  include  452,686 
tetrahedrons  (unstructured),  58,971  mixed  pyramids,  prisms  and  tetrahedrons  (transition), 
and  73,640  hexahedrons  (structured).  Note  that  while  the  unstructured  portion  is  1 1.4% 
of  the  tunnel  length,  it  accounts  for  87.4%  of  the  total  cells. 

Figure  3.16  shows  a  side  view  of  the  wing,  detailing  the  cell  face  distribution 
where  the  wing  interfaces  with  the  tunnel  wall.  Figure  3.17  displays  tightly  packed  cell 
faces  on  the  upper  wing  surface;  excluding  the  three  blowing  ports  and  the  region 
immediately  surrounding  them,  the  sides  of  each  triangular  face  measure  about  3.4  mm. 


Figure  3.16  -  Mesh  Generation,  Revision  L:  Profile  View  of  Wing  Mounted  to  Wall 
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Figure  3.18  is  a  close-up  view  of  one  of  the  blowing  ports  (enclosed  by  a  blue  circle)  on 

2 

the  wing  upper  surface,  where  104  triangular  faces  fit  into  the  0.49-mm  area. 


Figure  3.17  -  Mesh  Generation,  Revision  L:  Upper  Surface  of  Wing 


Figure  3.18  -  Mesh  Generation,  Revision  L:  Close-Up  View  of  Wing  Blowing  Port 
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Boundary  layer  resolution,  in  both  structured  and  unstructured  blocks  of  Revision 
L,  proved  adequate  for  some  turbulence  models  but  not  all.  Over  the  wing  surface, 
Figure  3.19  shows  that  wall  y+  values  ranged  between  15  and  60,  which  is  acceptable 
resolution  for  S-A,  k-s,  and  RSM  turbulence  models,  but  it  is  inadequate  for  the  LES 
model;  cells  around  the  blowing  ports  also  have  finer  resolution  as  shown  in  the  figure  by 
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Figure  3.19  -  Wall  v+  Values  along  Wing  Surface  in  x-Direction,  Revision  L,  Steady  S-A 
points  extending  as  low  as  v+  =  3.  Figure  3.20  shows  wall  v+  values  for  the  tunnel  walls. 
Structured  portions  of  the  domain  wall  boundary  layer  have  better  cell  resolution,  where 
v+  ranges  between  10  and  40,  and  unstructured  walls  have  poorer  resolution,  particularly 
in  the  transition  zones  between  structured  and  unstructured  regions  where  v+  reaches  as 
high  as  220  but  ranges  mostly  between  40  and  175.  This  indicates  that  the  unstructured 
tunnel  wall  resolution  is  acceptable  though  not  optimal;  finer  resolution  would  require  a 
tradeoff  for  increased  computing  time  and/or  resources.  Tunnel  wall  resolution  was 
generally  not  refined  in  this  study  because  emphasis  was  placed  on  solution  prediction  in 
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regions  of  vortex  generation  and  breakdown,  which  are  away  from  the  tunnel  walls; 
however  Revision  O  includes  some  boundary  layer  refinement  along  walls  and  wing. 
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Figure  3.20  -  Wally+  Values  along  Tunnel  Walls  in x-Direction,  Revision  L,  Steady  S-A 
As  for  additional  mesh  revisions,  Revision  M  is  specific  to  the  blowing  cases 
(Chapter  IV)  and  Revision  N  will  be  addressed  in  the  mesh  adaptation  section  later  in  this 
chapter.  Revision  O  addressed  the  boundary  layer  refresher  plate  approximations 
discussed  above.  Shown  in  Figure  3.21,  this  mesh  revision  altered  Revision  L  by 
creating  the  1.27-cm  gap  between  refresher  plate  and  tunnel  wall  and  did  not  include  the 
mounting  bracket,  seen  in  Figure  2.3;  the  plate  extends  to  the  top  and  bottom  tunnel  walls 
and  extends  downstream  45.72  cm  from  its  leading  edge,  as  in  the  experimental  setup. 
The  wing  remained  at  a  =  18°,  but  nodes  were  added  to  its  upper  surface  which  resulted 
in  an  additional  3,306  faces  and  58,730  tetrahedral  cells.  Nodes  were  also  added  to  the 
upstream  structured  block  to  refine  the  cell  transition  to  the  unstructured  volume  of  the 
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Figure  3.21  -  Mesh  Generation,  Revision  O:  Boundary  Layer  Plate  Separated  from 


Tunnel  Wall,  Plus  Wing  and  Transition  Refinements 
tunnel’s  wing  region  and  to  improve  y+  values  in  that  transition  zone;  this  resulted  in  an 
additional  14,188  hexahedral  and  14,731  mixed  cells.  Mesh  Revision  O  had  a  total  of 
207,784  nodes  and  672,946  volumetric  cells. 

These  modifications  did  improve  wall  y+  values  in  several  areas  within  the 
numerical  domain.  Lower  bounds  were  not  enhanced,  but  y+  upper  bounds  reduced  on 
the  wing  from  Revision  L’s  60  down  to  45;  along  the  tunnel  walls,  y+  reduced  its 
majority  upper  bound  (below  which  are  most  of  the  data)  from  175  down  to  100  and 
maximum  upper  value  from  220  down  to  130.  While  these  were  marked  improvements 
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in  cell  resolution  along  the  walls,  the  numerical  mesh  still  would  not  be  fine  enough  for 
the  LES  turbulence  model.  Further  results  from  this  grid  revision  are  included  later  in 
this  chapter,  but  it  did  not  prove  superior  to  mesh  Revision  L  because  approximating  the 
plate  as  flush  against  the  tunnel  wall  was  a  good  assumption. 

Solver  Parameters  and  Turbulence  Models 

This  section  discusses  the  procedure  used  in  FLUENT  to  initialize  the  various 
flow  cases,  including  selection  of  and  justification  for  solver  parameters,  flow  and 
turbulence  models,  discretization  schemes,  and  so  forth. 

Import  Numerical  Mesh 

After  opening  FLUENT  in  three-dimensional,  double-precision  mode,  the  user 
imports  or  reads  a  case  file,  which  in  this  instance  was  created  with  and  exported  by 
Gridgen  (using  Export  Analysis  Data  command).  A  Gridgen  case  file  contains 
information  about  the  mesh  and  boundary  conditions,  which  may  be  compatible  with 
FLUENT.  The  grid  was  then  checked  for  negative  volumes,  face  handedness  (where  left- 
handedness  on  a  face  indicates  negative  volume),  and  properly  matched  numbers  of 
nodes  to  cells  (FLUENT,  2001:  5.5.1).  Next  the  domain  was  reordered  and  cell  faces 
were  smoothed  and  swapped.  Reordering  the  domain  -  nodes,  faces  and  cells  -  places 
neighboring  cells  nearer  each  other  in  memory  “to  reduce  the  cost  of  memory  access” 
(FLUENT,  2001:  5.7.10),  and  smoothing/swapping,  which  applies  only  for  unstructured 
portions  of  a  grid,  repositions  nodes  and  combines  faces  where  possible  to  reduce  total 
cell  count  and  generally  improve  mesh  quality  (FLUENT,  2001:  23.11-23.11.2).  These 
actions  do  not  affect  the  model  geometry  proper.  Lastly,  the  case  file  in  FLUENT  was 
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scaled  from  its  default  meter  setting  to  the  centimeter  scale  used  in  generating  the  grid. 
All  other  units  remained  in  default  SI  units. 

Solver  Initialization 

Following  are  the  options  selected  in  FLUENT  to  generate  an  initial  solution.  In 
the  define^ models  menu,  solver  selections  included  coupled  solver,  implicit  fonnulation, 
three-dimensional  space,  and  steady  time.  The  coupled,  implicit  formulation  was  chosen 
for  computational  efficiency,  quicker  convergence  and  because  sufficient  computing 
resources  were  available.  Even  though  the  vortex  breakdown  problem  is  inherently 
unsteady,  a  steady  solver  was  chosen  initially,  hoping  that  it  would  lead  to  a  good  starting 
place  (mesh  geometry  determination  and  turbulence  model  selection)  before  running 
unsteady  cases.  For  later  cases  with  unsteady  time,  the  second-order  implicit  formulation 
was  selected;  second-order  refers  to  the  temporal  order  of  accuracy,  and  implicit  refers  to 
time-stepping  where  the  solution  iterates  to  convergence  at  each  time  step.  The  user  may 
opt  to  declare  a  maximum  number  of  iterations  per  time  step  vice  letting  it  run  to  meet  a 
set  of  convergence  criteria;  for  unsteady  cases  in  this  effort,  that  number  was  set  at  20, 
and  the  initial  time  step  was  set  at  0.0001  sec.  This  time  step  was  set  to  not  exclude 
effects  of  the  vortex  inner  core’s  observed  frequency  of  1,000  Hz  immediately  prior  to 
breakdown  (Novak,  Sarpkaya,  1999:  825);  it  was  later  discovered  that  this  time  step 
could  be  increased  to  0.0004  sec  without  adverse  effects  on  the  solution.  This  step  size, 
At  =  0.0004  sec,  allowed  essentially  for  2.5  computations  per  revolution  of  the  vortex 
core  upstream  of  breakdown.  The  solution  then  proceeded  through  an  undetermined 
number  of  time  steps,  but  at  least  a  few  flow  cycles,  until  the  solution  either  converged  or 
reached  a  relatively  steady  time-periodic  cycle.  One  flow  cycle  is  defined  as  the  time 
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required  for  an  average  particle  trace  to  flow  the  distance  of  the  domain,  or  domain  length 
divided  by  freestream  velocity;  for  Vx  =  16.05  m/s,  one  flow  cycle  is  0.436  sec. 

In  the  define— >  models  menu,  energy ;  equation  was  enabled  and  viscous  options 
include  inviscid,  laminar,  Spalart-Allmaras  (S-A),  k-s,  k-co  (not  investigated  in  this  study 
-  refer  to  Chapter  I),  Reynolds  Stress  (RSM),  and  Large-Eddy  Simulation  (LES).  Figure 
3.22  shows  the  FLUENT  define— >models^>viscous  drop-down  menu.  Initially,  the  S-A 
viscous  model  was  chosen  because  it  is  the  simplest  of  the  turbulence  models  available  in 
FLUENT  and  because  previous  studies  have  shown  that  turbulence  modeling  gives  more 


accuracy  in  prediction  of  vortex  strength,  bursting  location  and  extent  of  breakdown 
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region,  and  it  is  the  only  way  to  model  swirling  flow  with  anisotropic  turbulence 
(Ekaterinaris,  Schiff,  1990:  61;  Novak,  Sarpkaya,  1999:  833).  FLUENT  default  values 
were  used  for  S-A  model  coefficients,  and  no  attempts  were  made  to  alter  them.  Inviscid, 
laminar,  RNG  k-s  and  RSM  models,  with  FLUENT  default  settings,  were  all  compared  in 
steady  state  with  the  S-A  model.  The  following  models  were  run  in  unsteady  time:  LES, 
with  Smagorinsky-Lilly  subgrid  scale  model,  RNG  k-s,  RSM,  and  S-A  (FLUENT,  2001: 
10.4.2,  10.7.2).  With  exceptions  noted  above,  these  models  all  used  FLUENT  default 
settings  for  coefficients  and  wall  treatment. 

In  the  define^materials  menu,  the  only  operating  fluid,  air,  was  defined  with 
ideal-gas  density,  to  allow  for  compressibility  effects,  and  with  the  Sutherland  Law 
viscosity  using  Three-Coefficient  Method,  which  is  recommended  for  cases  with 
compressible  flow;  otherwise,  FLUENT  default  values  were  used  for  air  properties. 
Although  it  is  accurate,  Sutherland  viscosity  was  likely  overkill,  since  it  is  a  function  of 
the  ratio  of  total  and  static  temperature,  which  in  this  case  is  near  unity  (FLUENT,  2001: 
7.3.2).  Arguably  the  baseline  model  might  do  better  with  incompressible  air,  but  the  flow 
control  cases  inject  air  at  or  near  sonic  velocity,  so  the  more  robust  model  accounts  for 
compressibility  effects.  Under  define— ^operating  conditions,  the  operating  pressure  was 
set  to  zero  and  referenced  at  the  origin,  allowing  for  absolute  pressure  computations. 

Regarding  define— ^boundary  conditions,  the  wing,  unused  blowing  ports  on  the 
upper  wing  surface,  boundary  layer  refresher  plate,  and  all  tunnel  walls  were  defined  as 
no-slip,  aluminum,  impermeable,  insulated  wall  boundaries  (i.e.,  no  heat  or  mass  flux). 
In  the  far  field  mesh  models,  through  Revision  G,  remaining  domain  boundaries  were  set 
as  farfield  pressure,  where  gauge  pressure  was  99.56  kPa  (atmospheric  pressure  measured 
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in  the  lab  during  wind  tunnel  testing),  M*  was  0.0437  (based  on  the  original  claim  that  Vaj 
=  15  m/s),  temperature  was  298  K  (ambient  temperature  measured  during  tunnel  testing), 
turbulent  viscosity  ratio  (TVR)  was  10  -  FLUENT’S  recommendation  for  external 
compressible  flow  (Tutorial  3),  and  components  of  the  flow  direction  were  x  =  0.965926 
or  cos  15°,  y  =  0,  and  z  =  -0.258819  or  -sin  15°,  simulating  the  wing’s  angle  of  incidence. 
TVR  is  the  ratio  of  turbulent  to  laminar  viscosity,  where  a  value  of  unity  corresponds  to 
laminar  flow.  In  later  revision  cases,  it  was  discovered  that  varying  TVR  between  5  and 
10  had  no  significant  effect  on  the  solution  and  that  TVR  =  1  (in  farfield)  didn’t  allow  for 
sufficient  vortex  strength  (FLUENT,  2001:  27.4).  For  wind  tunnel  cases,  entry  to  the 
tunnel  was  defined  as  a  pressure  inlet,  where  initial  gauge  pressure  was  99.56  kPa,  total 
temperature  was  approximated  at  298  K,  TVR  was  unity  to  give  laminar  flow  at  the  inlet, 
flow  direction  was  normal  to  the  boundary,  and  gauge  total  pressure  was  set  to  a  value 
that  would  create  the  appropriate  freestream  velocity;  values  used  are  shown  as  Ptot  in 
Table  3.1.  Tunnel  exit  was  defined  as  a  pressure  outlet,  where  outlet  gauge  pressure  was 
99.56  kPa,  TVR  was  10  for  fully  turbulent  flow,  and  backflow  total  temperature  was  298 
K. 


Table  3.1-  Pressure  Inlet  and  Reference  Values  for  Various  Freestream  Velocities 


Voo  (m/s) 

Ptot  (Pa) 

P/v/t  (Pa) 

Poo  (kg/m3) 

Pdyn  (Pa) 

15.05 

99,692 

99,560 

1.164377 

131.867 

£ 

15.4 

99,698 

99,560* 

1.20767* 

143.206* 

16.05 

99,710 

99,560 

1.164437 

149.981 

17.03 

99,729 

99,560 

1.164500 

168.865 

19.99 

99,793 

99,560 

1.164713 

232.710 

24.98 

99,924 

99,560 

1.165149 

363.527 

T 


Conditions  reported  during  UC  wind  tunnel  testing 
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The  above  boundary  conditions  remained  the  same  whenever  possible  for  k-s, 
RSM  and  LES  turbulence  models;  for  inviscid  and  laminar  cases,  there  was  no  TVR  and 
Sutherland  viscosity  was  not  used  as  a  property  of  air.  For  k-s,  RSM  and  LES  models, 
turbulence  intensity  (ratio  of  root-mean-square  turbulent  fluctuations  magnitude  to 
freestream  velocity)  was  used  either  instead  of  or  in  addition  to  TVR,  where  intensity  was 
0.05%  at  the  pressure  inlet  and  10%  at  the  pressure  outlet  -  FLUENT’S  recommendations 
for  a  laminar  inlet  and  fully  turbulent  outlet  (27.4;  6.2.2). 

In  the  solve^controls  menu,  solution  parameters  were  set  for  discretization  type, 
Courant  or  CFL  number,  and  under-relaxation  factors.  Discretization  was  either  first-  or 
second-order  upwind  for  spatial  order  of  accuracy  and  for  time  accurate  cases  was  always 
second-order  implicit  for  temporal  order  of  accuracy.  CFL  numbers  ranged  from  5  to  10, 
based  on  whether  the  solution  was  stably  converging,  where  a  lower  CFL  number  gives  a 
more  stable  but  slower  converging  solution  and  a  higher  CFL  number  gives  quicker 
convergence  but  may  not  be  numerically  stable.  Under-relaxation  factors  remained  at  the 
FLUENT  default  settings  with  one  exception,  since  convergence  and  stability  problems 
were  typically  corrected  by  adjusting  spatial  order  of  accuracy  and/or  CFL  number;  an 
exception  was  with  the  RSM,  where  under-relaxation  factor  for  Reynolds  stresses  was 
changed  from  0.5  to  0.8  because  the  Reynolds  Stress  equation  variables  were  observed  to 
converge  quickly  and  stably. 

Once  models,  schemes,  discretization  methods,  boundary  conditions  and  so  forth 
were  judiciously  selected,  reference  values  were  established  (in  report  menu)  based  on 
pressure  farfield  or  inlet  conditions.  FLUENT-computed  reference  (or  freestream)  values 
included  density,  enthalpy,  pressure,  temperature,  velocity,  viscosity,  and  ratio  of  specific 
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heats;  pressure  (Pi nit  or  Poo),  velocity  and  density  values  that  were  used  are  indicated  in 
Table  3.1.  A  reference  value  which  required  user  specification  was  area;  for  this  case, 
that  area  was  one-half  the  planfonn  area,  0.033957  m  or  339.57  cm".  These  reference 
dimensions  and  values  were  then  used  with  pressure  values  on  the  wing  surface  to 
determine  the  integrated  quantities,  coefficients  of  lift  and  drag.  Next  the  domain  was 
partitioned  along  principal  axes  for  parallel  processing  on  between  6  and  20  processors; 
cases  with  larger  memory  and  computing  requirements,  typically  those  with  unsteady 
computations,  more  transport  equations  and/or  finer  mesh  resolution,  were  assigned  the 
greater  number  of  processors.  At  this  point,  the  domain  was  initialized  from  the  farfield 
or  tunnel  inlet  and  computational  iteration  commenced. 

Regarding  determination  of  how  many  partitions  to  create  within  the  model  for 
parallel  processing,  in  many  cases  a  tradeoff  was  made  between  reducing  computing  time 
and  having  enough  processors  available  to  allow  several  different  cases  to  run  in  concert 
with  each  other.  Another  consideration  was  whether  the  case  would  run  on  the  Aspen 
cluster  or  the  Hydra  cluster,  described  in  Chapter  II;  since  the  Aspen’s  processors  were 
slower  but  more  were  available,  its  cases  were  run  on  a  greater  number  of  processors  to 
give  roughly  equivalent  computing  time  per  iteration.  In  general,  steady  cases  ran  on  6,  8 
or  10  processors  and  were  partitioned  accordingly,  though  steady  cases  with  mesh 
adaptation,  and  hence  more  computational  cells,  ran  on  10  or  20  processors  -  10 
processors  when  other  cases  needed  to  run  and  20  when  resources  were  more  available. 
Steady  cases  converged  on  average  after  computing  for  8  to  12  wall-clock  hours;  mesh- 
adapted  steady  cases  required  between  15  and  38  hours  to  converge,  depending  on  the 
number  of  processors  and  number  of  cells  in  the  domain.  Unsteady  cases  ran  on  10 
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processors  when  other  cases  were  running  but  were  later  expanded  to  16  or  18  to 
accelerate  arrival  at  convergence.  Unsteady  cases  achieved  fully  developed  flow  after 
10-20  days  of  wall-clock  time.  As  for  a  limitation  on  number  of  partitions  in  the  domain, 
no  specific  study  was  performed. 

Computational  Plan  of  Attack 

In  general,  a  solution  with  second-order  spatial  and  temporal  accuracy  is  desired, 
and  this  was  attained  in  all  cases.  In  later  cases  the  solver  was  initialized  at  second-order, 
with  CFL  =  10  for  steady  cases  and  5  for  unsteady  cases,  and  successfully  generated  a 
second-order  solution.  Earlier  farficld  cases  required  running  the  first-order  solution  to 
convergence  then  changing  the  solver  to  second-order  and  running  that  solution  to 
convergence,  with  CFL  =  10  throughout.  Exceptions  included  the  k-e  and  RSM  steady 
cases  which  required  CFL  =  5  for  stability.  These  CFL  values  were  obtained  through 
trial  and  error  but  were  based  initially  on  FLUENT  recommendations  (ELUENT,  2001: 
Tutorials  3  and  4). 

Has  it  Converged  Yet? 

Convergence  is  attained  when  the  solution  no  longer  changes  with  more  iterations 
or  has  reached  a  fully  developed  flow  for  an  unsteady  case.  An  indicator  for  steady-state 
convergence  is  the  error  residual,  where  error  residuals  for  FLUENT’S  coupled  solver  are 
defined  as  the  square  root  of  the  averaged  time  rate  of  change  for  each  flow  variable 
within  the  domain  -  basically  the  difference  in  each  variable  from  one  iteration  to  the 
next.  Residuals  may  be  monitored  either  automatically  by  ELUENT  or  manually  by  the 
user  to  detennine  when  to  stop  iterating  the  solution.  A  general  guideline  for  steady 
convergence  is  when  the  residuals  have  decreased  by  three  orders  of  magnitude; 
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however,  this  does  not  necessarily  hold  true  for  turbulent  unsteady  cases,  where  it  was 
observed  that  some  residuals  reduced  by  as  many  as  seven  orders  of  magnitude  while 
others  reduced  only  by  three  (FLUENT,  2001:  22.16.1,  22.19.1). 

Occurrence  of  non-converging  residuals  necessitates  other  criteria  for 
convergence,  including  monitoring  for  steady  Cl,  Cd  and  outlet  mass  flow  rate  (or 
percent  change  in  mass  flow  rate).  Once  CL  and  Cd  have  little  or  no  change  from  one 
iteration  (steady)  or  time  step  (unsteady)  to  the  next,  the  solution  has  converged.  Also 
mass  flow  rate  may  be  monitored  at  the  outlet  surface  or  mass  flux  reports  may  be 
computed  to  determine  whether  mass  flow  is  conserved,  where  acceptable  mass  flow 
imbalance  should  be  no  more  than  0.5%  through  the  domain  (FLUENT,  2001:  22.16.3-4, 
Tutorial  4).  See  Figure  3.23  for  an  example  of  convergence  shown  by  steady  Cl  and  Cd. 


Drag  &  Lift  Convergence 


Jan  13.  20D3 

FLUENT  6.D  (3d.  dp,  coupled  imp,  S-A) 


Figure  3.23  -  Drag  and  Lift  Convergence  for  Revision  L,  a  =  18°,  Vm  =  16  m/s,  Spalart- 

Allmaras  Steady  Case 
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Another  possible  value  to  monitor  for  convergence,  though  at  one  specific  point, 


is  the  pressure  coefficient,  Cp.  It  is  calculated  as 


CP  = 


P~pi 


INIT 


(3.2) 


DYN 


where  p  is  local  static  pressure,  Pinit  is  reference  or  ambient  pressure,  and  Pdyn  is 
reference  dynamic  pressure. 


PDr„=  0.5 -p^-Vl  (3.3) 

with  tabulated  values  for  the  various  freestream  velocities  in  Table  3.1  (FLUENT,  2001: 
27.4).  In  this  study  Cp  was  not  monitored  for  solution  convergence,  but  its  distributions 
were  used  to  compare  solutions  from  different  meshes  and  turbulence  models  with  wind 
tunnel  test  data. 

For  all  steady  cases,  second-order  solution  convergence  was  achieved  with 
between  1,100  and  2,100  iterations.  For  unsteady  cases,  second-order  solution 
convergence  was  less  obvious.  In  these  cases,  the  error  residuals  reduced  overall  by 
between  three  and  seven  orders  of  magnitude  -  typically  one  to  two  orders  within  each 
implicit  time  step,  mass  flow  rate  was  conserved  to  within  a  few  thousandths  of  a  percent, 
but  while  Cl  and  C'd  values  developed  a  stationary,  periodic  cycle  in  time,  they  did  not 
reach  “steady”  values.  After  two  seconds  in  dimensional  time,  or  about  4.5  flow  cycles, 
the  flow  was  declared  fully  developed  since  average  Cl  and  C'd  values  were  then 
degrading/decreasing  by  only  0.4%  per  flow  cycle.  This  compared  to  acceptable 
variation  computed  in  a  similar  study,  where  CL  varied  by  1.6%  within  a  flow  cycle 
(Ekaterinaris,  Schiff,  1990:  61,  65).  At  a  step  size  of  0.0004  sec,  this  fully  developed  and 
mass-conserved  flow  required  5,000  time  steps  at  the  FLUENT-recommended  20 
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iterations  per  step  (FLUENT,  2001:  22.15.1),  for  a  total  of  100,000  iterations.  Within 
each  time  step,  the  20  iterations  generally  resulted  in  one  to  three  orders  of  magnitude 
reduction  in  error  residuals  for  the  flow,  continuity,  energy  and  transport  variables. 
Depending  on  the  number  of  processors  used  in  parallel,  fully  developed  unsteady 
solutions  required  between  278  (1 1.6  days)  and  472  (19.7  days)  wall-clock  hours. 

Determination  of  “Correct”  Initial  Conditions 

Studies  were  perfonned  to  determine  initial  conditions  which  would  give  an 
acceptable  solution  prediction.  Impetus  for  these  sub-studies  arose  from  uncertainty  in 
the  experimental  setup  and  resulting  data,  where  no  tolerance  was  initially  given  for  the 
freestream  velocity  measurement,  U,  =  1 5  m/s,  where  the  wing’s  angle  of  attack  was  not 
precisely  measured,  and  where  consideration  must  be  given  for  scraped  and  gouged  wing 
surfaces  and  a  leading  edge  subjected  to  light  grinding,  which  diminished  its  sharpness 
and  gave  it  a  slight  curvature  not  modeled.  As  mentioned  in  Chapter  I,  V was 
determined  to  be  15.4  m/s,  based  on  atmospheric  data,  but  without  benefit  of  a  tolerance 
or  error  bar.  These  factors  led  to  the  decision  to  adjust  the  wing’s  angle  of  attack  and  the 
freestream  velocity,  based  on  engineering  judgment,  to  compensate  for  testing 
uncertainty  and  to  give  a  closer  correlation  between  numerical  prediction  and 
experimental  data. 

For  these  evaluations  primary  consideration  was  given  to  how  closely  Cp  data 
compared  with  the  experiment  and  how  closely  Cl  compared  with  Polhamus’  theoretical 
value,  since  lift  had  not  been  experimentally  determined  for  this  setup.  Using  Equations 
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1.1,  1.3  and  1.4,  and  determining  Kp  and  Ky  coefficients  from  Polhamus’  plots,  for  a  = 
15°,  Cl  ~  0.7866  (Polhamus,  1971;  Polhamus,  1966). 

Selection  of  Freestream  Velocity 

Results  from  earlier  cases,  using  mesh  Revisions  H  and  K,  showed  that  varying 
freestream  velocity  (increasing  it  from  15  m/s  to  17  or  25  m/s)  did  not  affect  Cp 
predictions  because  it  increased  magnitude  of  dynamic  and  local  pressures  in  proportion 
to  each  other;  increased  velocity  also  did  not  significantly  change  Cl,  where  the  increase 
to  25  m/s  resulted  in  a  1.7%  increase  in  lift  coefficient.  Thus  Vx=  16  m/s  was  chosen 
because  it  increased  local  static  pressures  closer  to  the  experimentally  detennined  values 
without  significantly  increasing  the  reference  dynamic  pressure  (150  Pa)  above  the 
average  measured  value  (143  Pa)  and  because  it  was  within  a  reasonable  tolerance  of  the 
estimated  15.4  m/s  wind  tunnel  condition. 

Selection  of  Wing  Angle  of  Attack 

Using  mesh  Revision  L,  cases  were  run  for  angles  of  attack,  a  =  15,  16,  18,  19,  20 
and  21°,  using  the  Spalart-Allmaras  (S-A)  turbulence  model  in  steady  state  and  Vx  = 
16.0467  m/s.  Cl  values  from  several  experiments,  from  this  study,  and  from  Polhamus 
are  listed  in  Table  3.2.  The  “Experiment”  Cl  value  came  from  testing  performed  at 
Louisiana  State  University  several  years  ago,  where  the  same  wing  was  attached  to  a 
half- fuselage.  Those  results  showed  that  such  a  configuration  creates  an  effective 
increased  angle  of  attack  so  Cl  for  the  wing  without  fuselage  would  be  somewhat  less, 
likely  less  than  the  theoretical  value  (Guillot,  1999:  33).  In  Chapter  I  and  as  shown  by 
empirical  data  from  Wentz,  Kohlman,  Seginer  and  Salomon  in  Table  3.2,  it  was 
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established  that  Polhamus’  estimation  tended  to  be  high,  so  cases  were  eliminated  with 
Cl  greater  than  about  0.79,  leaving  as  options  cases  for  a  =  15,  16  and  18°.  Note  also  that 
numerically  estimated  Cl  values  correspond  closely  with  experimental  data  in  the  table; 
this  validates  the  steady,  S-A  numerical  model’s  ability  to  accurately  predict  Cl- 

Table  3.2  -  Comparison  of  Lift  Coefficient  Resulting  from  Various  Angles  of  Attack 


Case 

a  (deg) 

Vao  (m/s) 

Experiment a 

~15 

-15.4 

0.81 

Theory  (Polhamus) 

15 

~0 

0.7866 

18 

~0 

0.9614 

Wentz,  Kohlman b 

15 

18 

0.72 

18 

18 

0.80 

Seginer,  Salomon c 

15 

30 

0.70 

18 

30 

0.84 

Revision  L  (S-A) 

15 

16.0467 

0.6851 

16 

66 

0.7240 

18 

66 

0.7974 

19 

66 

0.8268 

20 

66 

0.8587 

21 

66 

0.8812 

3  Results  with  fuselage  attached  to  wing  (Guillot,  1999:  32) 

Delta  wing  with  A=60°,  sharp  leading  edge  but  only  20%  of  wing  thickness 


in  current  study  (Wentz,  Kohlman,  1971:  157-158) 

Delta  wing  with  A=60°  and  relatively  sharp  leading  edge  (Seginer, 
Salomon,  1986:  803-804) 


Other  factors  affecting  lift  coefficient  include  whether  the  wing  was  mounted 
flush  against  the  tunnel  wall,  whether  inlet  flow  was  laminar  or  turbulent,  and  whether 
boundary  layer  refreshing  was  done.  The  difference  between  mesh  Revisions  J  and  K, 
for  same  angle  of  attack,  was  that  the  wing  was  flush  against  the  wall  in  Revision  K;  this 
increased  Cl  by  5.0%.  With  Revision  K,  two  cases  were  run  with  TVR  set  at  1  and  10, 
ah  other  conditions  being  the  same;  the  case  with  TVR  =  10,  or  turbulent  inlet  flow, 
increased  Cl  by  0.6%.  Mesh  Revision  J  differed  from  Revision  H  in  that  it  simulated 
boundary  layer  refreshing  by  removing  nearly  2  m  of  tunnel  wall  upstream  of  the  wing; 
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this  change  resulted  in  increasing  Cl  by  10.1%.  Neglecting  the  insignificant  difference 
created  by  varying  inlet  turbulence,  Revision  L  incorporated  both  mounting  wing  flush 
with  the  wall  and  introducing  a  boundary  layer  refresher  plate. 

The  observations  above  did  not  conclusively  indicate  which  angle  of  attack 
should  be  used  in  the  numerical  model,  since  Cl  data  were  not  available  from  testing  in 
this  configuration.  Therefore,  the  following  detailed  comparisons  of  CP  data,  Figures 
3.24-27,  were  used  instead,  which  indicated  that  the  best  option  for  angle  of  attack  was  a 
=  18°. 

Figure  3.24  contains  data  along  the  chordwise  location,  x/c  =  0.35,  from  the 
experiment  conducted  at  University  of  Cincinnati,  seen  as  a  blue  line  with  a  square 
marker  for  each  pressure  sensor  location  in  the  spanwise  direction.  Each  experimental 
data  point  was  calculated  using  Equation  3.2,  where  p  was  pressure  measured  at  that 
specified  location,  Pinit  was  ambient  pressure  in  the  lab,  and  Pdyn  was  dynamic 
freestream  pressure  measured  in  the  wind  tunnel.  This  resulted  in  negative  pressure 
coefficients,  but  they  were  plotted  on  a  positive  scale  for  easier  viewing,  as  is  common 
practice.  Numerical  solutions,  at  the  various  angles  of  attack,  are  shown  as  lines  with  no 
markers  except  for  the  a  =  18°  case,  which  is  indicated  as  a  black  line  with  markers. 
Numerical  data  points  corresponding  to  predicted  static  pressure  values  were  extracted 
from  the  solution  along  the  spanwise  line  at  chordwise  location,  x/c  =  0.35,  and  were 
plugged  into  Equation  3.2  to  yield  the  according  Cp.  This  procedure  was  performed 
instead  of  using  the  FLUENT-computed  Cp  values  to  allow  greater  semblance  to  the 
procedure  for  calculating  experimental  CP.  Figures  3.25-27  provide  data  for  x/c  =  0.55, 
0.75  and  0.95,  respectively. 
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In  the  following  observations,  numerical  cases  with  a  =  15  and  16°  are  referred  to 
as  “lower-alpha  cases,”  and  a  =  19,  20  and  21°  are  “higher-alpha  cases.”  In  Figure  3.24, 
a  =  18°  is  the  best  case  because  its  slope  is  better  than  those  of  the  higher-alpha  cases  and 
its  magnitude  is  greater  than  the  lower-alpha  cases.  While  its  magnitude  is  lower  than  the 
higher-alpha  cases,  it  maintains  better  Cl  correlation.  The  same  arguments  hold  for  the 
remaining  figures,  and  note  in  Figures  3.26  and  3.27  that  the  a  =  18°  case  has  the  greatest 
peak  value  and  has  a  slope  which  closely  matches  that  of  the  experimental  data.  Of 
prime  importance  is  the  model’s  ability  to  closely  approximate  the  solution  in  the  region 
inside  and  following  vortex  breakdown  at  or  near  x/c  =  0.4,  and  the  a  =  18°  case  indeed 
appears  to  do  this.  Furthermore,  this  case  clearly  predicts  pressure  values  with  greater 
accuracy  than  the  a  =  15°  case. 


Figure  3.24  -  Angle  Selection  Data  at  x/c  =  0.35  (Mesh  Revision  L,  Steady  S-A) 
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Figure  3.25  -  Angle  Selection  Data  at  x/c  =  0.55  (Mesh  Revision  L,  Steady  S-A) 


Angle  Selection  -  Experimental  and  CFD  Model  (RevL,  S-Asteady) 
x/c  =  0.75 


Figure  3.26  -  Angle  Selection  Data  at  x/c  =  0.75  (Mesh  Revision  L,  Steady  S-A) 


3-38 


Angle  Selection  -  Experimental  and  CFD  Model  (RevL,  S-A steady) 
x/c  =  0.95 


r 


A 


\ _ J 

Figure  3.27  -  Angle  Selection  Data  at  x/c  =  0.95  (Mesh  Revision  L,  Steady  S-A) 

Therefore,  to  compensate  for  uncertainties  in  experimental  measurement  and 
procedure,  this  study’s  final  or  best  numerical  mesh  and  flow  model  set  the  wing  at  18° 
angle  of  attack  and  freestream  velocity  at  16  m/s.  It  must  be  pointed  out  that  these  sub¬ 
studies  were  performed  using  the  steady  solver  with  S-A  turbulence  model;  the  section  on 
time  accurate  solutions  later  in  this  chapter  shows,  however,  that  more  accurate  solutions 
(i.e.,  those  which  more  precisely  model  the  highly  unsteady  nature  of  the  flow)  used  an 
unsteady  fonnulation.  However,  since  the  steady  S-A  results  did  show  good  correlation 
with  experimental  findings,  the  initial  conditions  above  were  deemed  acceptable. 
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Results 


Effect  of  Model  Geometry 

The  initial  numerical  mesh  through  Mesh  Revision  F  included  the  wing,  offset 
from  a  no-slip  wall  boundary,  and  a  pressure  far li eld  which  extended  to  various 
dimensions  from  the  wing.  It  was  likely  a  flawed  approach  to  include  a  constraining  wall 
in  this  manner,  particularly  since  that  wall  eventually  exceeded  the  dimensions  of  the 
equivalent  adjacent  wall  in  the  wind  tunnel  experiment.  Better  options  could  include 
modeling  the  experimental  setup,  as  was  done  in  ensuing  mesh  revisions,  or  to  establish  a 
symmetry  plane  instead  of  a  wall  in  the  far  lie  Id  cases,  where  a  symmetry  plane  would 
ensure  no  wall  boundary  layer  influence  on  the  wing.  However,  a  case  with  symmetry 
plane  would  resemble  even  less  the  actual  experimental  setup,  thus  that  option  was  not 
exercised.  The  following  Cp  plots,  compared  with  experimental  data,  Figures  3.28-31, 
show  that  none  of  the  farficld  models  adequately  predicted  the  experimental  results. 

For  each  solution,  the  steady  case  used  the  S-A  turbulence  model,  a  =  15°  and  Vx 
=  15  m/s,  largely  for  the  sake  of  consistent  evaluation;  while  none  of  these  solutions 
could  be  precise  with  a  steady  solver  (refer  to  discussion  in  Chapter  I),  relative 
comparison  was  sufficient  for  establishing  the  most  suitable  model  with  which  further 
investigation  could  proceed.  As  an  additional  note  on  using  a  steady  solver  to  ascertain 
the  best  configuration  to  run  the  unsteady  cases:  steady  solutions  predicted  CP  values 
closer  to  experimental  values  than  did  any  unsteady  solution,  because  steady  or  time 
averaged  solutions  more  precisely  represent  the  experimental  data,  which  were  also  time 
averaged;  however,  to  reiterate,  steady  solvers  cannot  accurately  predict  all  the  flow 
physics  involved  with  vortex  generation  and  particularly  vortex  breakdown.  Thus  the 
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preferred  model  should  be  able  to  adequately  predict  Cp  data  and  physics  of  flow, 
specifically  vortex  strength  (velocity  and  pressure)  and  location  of  vortex  breakdown. 


Figure  3.28  -  Cp  Comparison  of  Initial  through  Mesh  Revision  F  at  x/c  =  0.35 
(Steady  S-A,  a  =  15°,  =  15  m/s) 


In  Figures  3.28-31,  Cp  predictions  generated  by  the  solution  from  the  initial 
numerical  mesh  are  indicated  by  a  red  line,  experimental  data  are  indicated  by  a  blue  line 
with  square  markers,  and  predictions  from  mesh  Revisions  A-F  are  indicated  by  green, 
blue  (no  markers),  gray,  black,  purple  and  light  blue,  respectively.  In  each  of  the  plots, 
the  peak  values  (actually  a  minimum  negative  value)  indicate  the  experimentally 
determined  or  numerically  predicted  center  location  of  the  vortex  core,  since  a  large 
negative  coefficient  represents  a  static  pressure  lower  than  ambient  conditions,  which 
corresponds  also  to  an  area  of  higher  velocity.  The  initial  mesh  tended  to  result  in 
overestimated  Cp  predictions  but  was  reasonably  close  in  prediction  of  vortex  peak 
location,  because  the  initial  numerical  model  gave  no  offset  between  wing  and  wall. 


3-41 


Revisions  A-F  predicted  vortex  peak  locations  which  strayed  from  the  experimentally 
determined  location  by  3-6%.  This  likely  corresponds  to  the  numerical  wing  being  offset 
from  the  wall  by  1.27  cm,  which  is  6.1%  of  the  half  span,  representing  an  effectual 
decrease  in  angle  of  attack  or  an  increase  in  sweep  angle.  In  Figure  3.28,  the  early 
Revisions  -  B,  C  and  D  -  provided  the  greatest  relative  improvement  by  enlarging  the 
numerical  domain  downstream  and  spanwise,  while  Revisions  E  and  F  resulted  in  worse 
predictions  because  their  upstream  expansions  resulted  in  greater  boundary  layer 
generation  along  the  wall,  thus  negating  the  experiment’s  boundary  layer  refreshing. 
Note  similar  trends  in  Figures  3.29-31.  At  x/c  =  0.35,  Revisions  A-F  predicted  peak 
values  with  an  error  between  34  and  49%  of  the  experimentally  detennined  peak,  which 
is  of  course  unacceptable.  These  Revisions’  errors  were  about  10-33%  at  x/c  =  0.55,  2- 
23%  at  x/c  =  0.75,  and  1-21%  at  x/c  =  0.95.  These  results  would  be  marginally 
acceptable  if  Cp  magnitude  were  the  only  criterion,  but  that  was  not  the  case. 

All  of  the  numerical  meshes  result  in  better  pressure  prediction  in  these  plots  as 
the  flow  advances  downstream  over  the  wing  surface;  this  is  likely  because  vortex  energy 
increases  with  increasing  turbulence,  and  as  a  result  pressure  decreases  and  gives  a  larger 
negative  pressure  coefficient.  It  may  also  be  a  function  of  numerical  grid  resolution, 
where  finer  resolution  may  be  required  near  the  apex  to  resolve  the  smaller  vortex 
diameter  and  higher  rotational  velocity.  While  it  appears  that  Revision  B  gave  the  closest 
estimations,  its  domain  did  not  contain  the  flow  solution;  also  while  Revisions  D,  E  and  F 
did  contain  the  solution  (extent  of  vortex)  within  the  computational  domain,  they  either 
failed  to  provide  domain  extents  acceptable  for  low  subsonic  flow  (as  with  Revision  D, 
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Figure  3.29  -  Cp  Comparison  of  Initial  through  Mesh  Revision  F  at  x/c  =  0.55 
(Steady  S-A,  a  =  1 5°,  F»  =  1 5  m/s) 


which  needed  more  upstream  extent)  or  failed  to  provide  results  within  an  acceptable 
margin  of  the  experimental  data.  Again,  none  of  Revisions  Initial  through  F  adequately 
predicted  vortex  location  or  Cp  values  at  x/c  =  0.35.  Farther  downstream  on  the  wing 
surface,  these  models  predicted  closer  Cp  magnitudes  because  flow  became  more 
turbulent  in  this  region,  but  they  still  failed  to  predict  the  correct  vortex  location.  This 
lack  of  improvement  in  predictions  led  to  creation  of  numerical  meshes  which  included 
the  wind  tunnel  (Revisions  G-J)  and  boundary  refresher  plate  (Revisions  K-O). 
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Figure  3.30  -  Cp  Comparison  of  Initial  through  Mesh  Revision  F  at  x/c  =  0.75 
(Steady  S-A,  a  =  15°,  ¥«>  =  15  m/s) 


Baseline  Comparison  -  Experimental  and  CFD  Initial  Model  thru  RevF 
x/c  =  0.95 


Figure  3.3 1  -  Cp  Comparison  of  Initial  through  Mesh  Revision  F  at  x/c  =  0.95 
(Steady  S-A,  a  =  15°,  =  15  m/s) 
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To  review  the  changes  implemented  in  the  next  mesh  revisions:  Revision  G 
modeled  the  wind  tunnel  with  wing  offset  from  the  wall  and  with  tunnel  inlet  2  m 
upstream  of  the  wing  (results  were  not  obtained  for  Revision  G,  since  it  was  revised 
before  a  converged  solution  was  obtained);  Revision  H  improved  resolution  in  the 
boundary  layer  region  and  maintained  the  wing’s  offset  from  the  wall  and  the  tunnel  inlet 
location;  Revision  I  maintained  the  wing  offset  but  truncated  the  tunnel  inlet  location  to 
about  two  chords  upstream  of  the  wing  apex  to  assess  the  impact  of  boundary  layer  flow 
over  the  wing;  Revision  J  still  maintained  the  wing  offset  from  the  tunnel  wall  and 
reduced  the  inlet  distance  to  several  centimeters  upstream  of  the  wing  apex,  the  first 
intentional  attempt  to  simulate  boundary  layer  refreshing;  Revision  K  placed  the  wing 
flush  against  the  tunnel  wall  and  corrected  the  distance  (shorter)  between  wing  apex  and 
leading  edge  of  refresher  plate;  Revision  L  included  the  refresher  plate’s  45°  ramp  or 
bevel  and  replaced  the  tunnel  inlet  to  2  m  upstream  of  the  wing  apex;  and  Revision  O 
modeled  the  wing  flush  against  the  refresher  plate,  which  plate  was  then  offset  from  the 
tunnel  wall,  and  tunnel  inlet  at  2  m  from  the  wing.  Revision  O  came  as  an  afterthought  to 
validate  or  refute  the  assumption  of  placing  the  refresher  plate  flush  against  the  tunnel 
wall,  and  its  evaluation  showed  the  assumption  for  Revision  L  was  good  and  allowed  for 
a  simpler  numerical  geometry. 

All  of  the  above  meshes  modeled  5  m  of  tunnel  downstream  of  the  wing,  and  in 
the  following  comparisons  used  the  steady  solver,  S-A  model,  a  =  15°  and  Vx  =  15  m/s, 
except  Revision  L  used  Vx=  16  m/s.  (As  noted  earlier  in  this  chapter,  the  slightly  higher 
freestream  velocity  showed  no  effect  on  Cp  values.)  Also  the  comparison  between  Mesh 
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Revisions  L  and  O  used  a  =  18°,  F»  =  16  m/s,  and  the  steady  RSM  turbulence  model, 
since  the  S-A  model  introduced  instability  in  the  solver  for  Revision  O. 

In  Figures  3.32-35,  data  at  the  four  chordwise  locations  (x/c  =  0.35,  0.55,  0.75  and 
0.95)  are  shown  by  a  blue  line  with  square  markers  for  the  experiment  and  by  lines  with 
no  markers  of  colors,  green,  blue,  gray,  black  and  red  for  predictions  from  Mesh 
Revisions  H-L,  respectively. 

In  these  plots,  Revisions  K  and  L  give  the  closest  prediction  of  peak  Cp  location, 
which  corresponds  to  vortex  core  centerline,  because  they  both  model  the  wing  as  flush 
against  the  wall;  Revisions  H-J  incorrectly  model  the  wing  as  offset  from  the  wall. 
Revision  L  predicts  closer  than  does  Revision  K  to  the  experimentally  determined 
centerline  at  all  chordwise  locations  except  x/c  =  0.35,  where  they  predict  the  same 
location.  This  is  because  Revision  L  resembles  more  closely  the  actual  experimental 
setup,  which  includes  boundary  layer  buildup  from  the  2  m  of  tunnel  wall  upstream  of  the 
wing. 

While  some  discussion  of  Cp  magnitude  follows,  primary  consideration  must  be 
given  to  the  model  which  most  closely  predicts  vortex  core  location  because  development 
of  this  model  is  intended  ultimately  to  be  part  of  a  control  system  which  needs  to  direct 
blowing  momentum  into  the  vortex  core.  To  do  this,  the  vortex  core  location  must  be 
accurately  predicted.  Since  these  data  are  all  time  averaged,  including  the  experiment, 
except  for  the  row  of  sensors  upstream  of  vortex  breakdown  (at  x/c  =  0.35)  none  of  them 
shows  a  true  location  of  the  vortex  core,  since  it  meanders  back  and  forth  in  time  in  a 
spanwise  direction  once  vortex  bursting  has  occurred  (see  Chapter  I  for  more  infonnation 
and  references);  this  phenomenon  is  shown  in  later  unsteady  cases.  Nonetheless,  steady 
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modeling  provided  relatively  accurate  predictions  from  the  varied  numerical  meshes  and 
thus  serves  to  identify  the  best  numerical  mesh. 

In  these  comparisons,  little  consideration  is  extended  to  actual  predicted  Cp 
magnitudes  because  they  were  subsequently  altered  to  more  closely  match  experimental 
data  by  adjusting  wing  angle  of  attack;  hence  relative  comparisons  are  made  to  detennine 
the  best  computational  grid.  Revision  H  clearly  predicted  the  least  accurate  solution  at  all 
locations.  Also  at  each  location,  Revisions  I  and  J  which  progressively  truncated  the 
distance  between  wing  apex  and  tunnel  inlet,  showed  marked  improvement  over  Revision 
H  results.  This  indicated  the  clear  need  for  boundary  layer  refreshing  in  the  numerical 
model.  Revision  K’s  results  showed  the  combined  effect  of  moving  the  wing  to  be  flush 
with  the  wall  and  again  reducing  distance  to  the  tunnel  inlet.  Thus  at  each  chordwise 
position,  the  Cp  peak  moved  closer  to  that  of  the  experimental  data,  and  Cp  magnitudes 
increased  at  x/c  =  0.35  and  0.55.  For  x/c  =  0.75  and  0.95,  Revision  K’s  predicted  Cp 
magnitudes  increased  over  those  of  Revision  J,  as  shown  moving  a  spanwise  outward 
direction  in  Figures  3.34  and  3.35,  but  then  it  predicts  lower  peak  values  around  the 
vortex  core.  This  indicates  the  advantage  of  a  fresher  boundary  layer  was  then 
overridden  by  boundary  layer  or  other  viscous  effects  from  having  the  wing  flush  against 
the  tunnel  wall  and  from  interaction  with  flow  over  the  other  tunnel  walls.  This  also 
showed  that  wall  or  boundary  layer  effects  from  the  tunnel  have  a  diminishing  influence 
on  prediction  of  vortex  strength  in  the  region  aft  of  vortex  breakdown.  Revision  K, 
however,  was  superior  to  Revision  J  because  it  better  predicted  vortex  centerline  location. 
Revision  L’s  introduction  of  the  boundary  layer  plate  and  extended  upstream  tunnel,  as 
opposed  to  approximating  by  a  closer  tunnel  inlet,  resulted  in  decreased  vortex  strength 
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but  more  accurate  vortex  core  position.  Clearly,  structure  of  flow  approaching  the  wing 
is  one  of  the  more  significant  factors  affecting  solution  accuracy.  Therefore,  Mesh 
Revision  L  was  selected  for  further  investigations  because  it  most  closely  predicted 
vortex  centerline  position  and  it  most  closely  modeled  the  experimental  setup. 


Figure  3.32  -  Cp  Comparison  of  Mesh  Revision  H  through  L  at  x/c  =  0.35 
(Steady  S-A,  a  =  15°,  F®  =  15  m/s  (except  for  RevL,  F®  =  16  m/s)) 


3-48 


Figure  3.33  -  CP  Comparison  of  Mesh  Revision  H  through  L  at  x/c  =  0.55 
(Steady  S-A,  a  =  15°,  Vx  =  15  m/s  (except  for  RevL,  Vm  =  16  m/s)) 


Figure  3.34  -  Cp  Comparison  of  Mesh  Revision  H  through  L  at  x/c  =  0.75 
(Steady  S-A,  a  =  15°,  Voo  =  15  m/s  (except  for  RevL,  Foo  =16  m/s)) 
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Figure  3.35  -  Cp  Comparison  of  Mesh  Revision  H  through  L  at  x/c  =  0.95 
(Steady  S-A,  a  =  15°,  Vx=  15  m/s  (except  for  RevL,  Vx  =  16  m/s)) 

To  lend  further  validity  to  the  assertion  that  Mesh  Revision  L  was  superior  to 


Revision  K,  Figure  3.36  compares  experimental  and  predicted  CP  values,  where  all  four 
chordwise  data  sets  are  captured  in  one  plot  and  where  solutions  from  Revisions  K  and  L 
were  computed  with  a  =  18°  and  Vx  =  16  m/s.  These  data  show  that  while  Revision  K 
does  better  at  x/c  =  0.35,  Revision  L  predicts  a  solution  comparable  or  better  at  the  other 
chordwise  locations. 
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Baseline  Compare  -  Exper  (blue  sq)  and  CFD  Model  RevK  (red)  and  Rev  L  (green),  Steady  S-A 
left  to  right  -  x/c  =  0.35,  0.55,  0.75,  0.95 


r 


,  y/s  (percent) 

Figure  3.36  -  CP  Comparison  of  Revisions  K  and  L  (Steady  S-A,  a  =  18°,  F»  =  16  m/s) 
A  final  consideration  in  determining  correct  numerical  model  geometry  for  the 
flow  model  was  to  evaluate  the  effect  of  the  boundary  layer  refreshing  plate’s  placement 
-  flush  with  or  offset  from  the  tunnel  wall.  Using  the  steady  RSM  turbulence  model, 
Revision  L’s  solution  computed  Cl  =  0.740  and  Revision  O  gave  Cl  =  0.768,  thus 
Revision  O  resulted  in  a  Cl  prediction  closer  to  theoretical  and  experimental  values, 
shown  earlier  in  Table  3.2.  Data  in  Figure  3.37  show  that  Mesh  Revision  O’s  (green  line) 
geometry  gives  slight  improvement  in  CP  magnitude  and  slope  over  that  of  Revision  L; 
however,  using  the  steady  RSM  solver,  neither  computational  mesh  resulted  in  close 
comparison  with  experimental  results.  In  addition  to  comparing  Cp  and  Cl  values,  a 
qualitative  comparison  between  boundary  layers  follows  in  Figures  3.38  and  3.39. 
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Baseline  Compare  -  Exper  (blue  sq)  and  CFD  Model  RevL  (red)  and  RevO  (green),  Steady  RSM 
left  to  right-  x/c  =  0.35,  0.55,  0.75,  0.95 


y/s  (percent) 

Figure  3.37  -  CP  Comparison  of  Revisions  L  and  O  (Steady  RSM,  a  =  18°,  Vx  =  16  m/s) 


In  Figures  3.38  and  3.39,  velocity  predictions  were  nondimensionalized  with  the 


freestream  velocity  then  truncated  at  V/V„  =  0.50,  where  V  is  velocity  magnitude,  for 


cleaner  viewing  of  the  layer.  The  only  noteworthy  difference  between  these  two  figures 


is  that  Figure  3.39  shows  formation  of  a  thicker  boundary  immediately  upstream  of  the 


refresher  plate  leading  edge,  which  thickness  then  disappears  into  the  gap;  otherwise,  the 


boundary  layer  flow  approaching  the  wing  in  each  case  is  nearly  identical  with  or  without 


an  offset  between  the  refresher  plate  and  tunnel  wall. 
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Figure  3.38  -  Boundary  Layer  Effects  with  Refresher  Plate  Flush  with  Tunnel  Wall 
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Figure  3.39  -  Boundary  Layer  Effects  with  Refresher  Plate  Offset  from  Tunnel  Wall 
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Figure  3.40  shows  via  streamlines  that  flow  around  the  boundary  layer  refresher 
plate  is  laminar  and  steady,  unless  or  until  the  flow  rolls  into  vortex  effects,  as  observed 
with  the  purple  streamlines.  The  black  streamlines  indicate  that  flow  over  the  refreshing 
plate  upper  surface  is  undisturbed  until  it  mixes  with  the  downstream  turbulence  from  the 
vortex  and  its  breakdown,  shown  with  blue  streamlines.  While  Mesh  Revision  O, 
compared  with  Revision  L,  apparently  gives  slightly  more  accurate  Cp  and  Cl  predictions 
and  does  not  overly  complicate  the  flow  field,  its  computation  required  15.4%  more  time 
per  iteration  and  20%  more  memory.  Therefore,  results  using  the  modeled 
approximations  of  Mesh  Revision  L  were  deemed  acceptable,  and  the  assumption  of 
placing  the  refresher  plate  flush  against  the  tunnel  wall  was  appropriate. 


Figure  3.40  -  Streamlines  Showing  Benign  Effect  of  Refresher  Plate  Offset  from  Wall 
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In  summary,  Mesh  Revision  L  was  selected  as  an  acceptable  numerical  mesh 
because  it  modeled  correct  tunnel  geometry,  it  made  a  suitable  approximation  by  placing 
the  refresher  plate  flush  against  the  tunnel  wall,  it  resulted  in  accurate  prediction  of 
vortex  core  position,  and  its  predicted  Cp  and  Cl  values  (with  modified  wing  angle  of 
attack)  were  within  a  reasonable  range  of  experimentally  detennined  values. 

Effect  of  Turbulence  Model  with  Time  Averaged  Flow 

While  time  averaged  or  steady  flow  solutions  may  not  accurately  predict  all  the 
physics  of  this  highly  unsteady  flow,  they  do  provide  cursory  correlation  with  the 
experimental  data,  where  each  Cp  point  on  the  plots  below  represents  the  average  of  500 
collected  data  points  at  each  respective  pressure  sensor.  Thus  time  averaged,  numerically 
predicted  Cp  data  allow  for  the  most  correct  correlation  between  results,  whereas  time 
accurate  solutions  show  snapshots  in  time  and  not  an  averaged  value. 

Figures  3.41-44  show  experimental  Cp  data  compared  with  predicted  values  using 
inviscid  and  laminar  (non-turbulence)  models,  and  S-A,  RNG  k-s,  and  RSM  turbulence 
models.  None  of  the  numerical  models  showed  good  agreement  with  the  experiment,  but 
the  steady  S-A  model  came  closest  to  predicting  a  correct  solution,  with  respect  to  Cp 
correlation.  It  accurately  predicted  the  vortex  core  centerline  location  at  all  chordwise 
positions,  and  it  came  closest  to  slope  and  peak  magnitude  at  x/c  =  0.55  and  0.75,  which 
are  both  aft  of  the  vortex  bursting  location  -  a  region  where  accurate  prediction  is  crucial. 
At  x/c  =  0.95,  the  S-A  model  best  predicted  the  slope  or  region  of  vortex  activity,  and  the 
magnitude  was  closer  than  the  other  turbulence  models. 

With  some  modification  to  help  it  predict  the  physics  of  vortex  bursting,  S-A  may 
prove  to  be  the  best  of  these  turbulence  models.  One  modification  (termed  S-A  Rotation 
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Correction)  allows  “production  of  turbulent  viscosity  to  be  reduced  in  regions  of  high 
vorticity,”  and  another  (termed  Detached  Eddy  Simulation)  combines  the  computational 
advantages  of  the  Reynolds-Averaged  Navier-Stokes  method  with  the  accuracy  of  LES  in 
computing  highly  separated  flow  fields  (Morton,  Forsythe,  Mitchell,  Hajek,  2002:  3,  5). 
But  for  this  study,  the  S-A  model  proper  remains  inadequate  since  it  failed  to  correctly 
predict  occurrence  of  vortex  breakdown. 

At  x/c  =  0.35,  all  of  the  numerical  models  predicted  nearly  the  same  solution, 
though  the  laminar  case  best  estimated  CP  slope  and  magnitude,  likely  because  this 
location  most  closely  corresponds  to  laminar  flow  -  prior  to  vortex  breakdown.  At  this 
location,  no  model  accurately  predicted  vortex  strength  or  size,  which  may  have  been  due 
to  coarse  grid  resolution,  shown  in  greater  detail  in  the  next  section.  The  laminar  model’s 
predicted  Cp  curves  lost  smoothness  after  vortex  breakdown  occurred,  seen  in  Figures 
3.42-44,  indicating  that  a  laminar  model  unsurprisingly  has  difficulty  predicting 
properties  of  highly  turbulent  flow.  The  inviscid  solution  coincided  relatively  closely 
with  laminar  predictions,  but  since  the  flow  field  was  highly  dependent  upon  freestream 
conditions  over  the  wing  and  since  those  conditions  included  no  boundary  layer  flow,  the 
inviscid  prediction  is  not  a  valid  option  for  further  investigation.  It  was  included  merely 
for  pedagogical  rigor. 
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Figure  3.41  -  CP  Comparison  of  Steady  Turbulence  Models  at  x/c  =  0.35 
(Mesh  Revision  L,  a  =  18°,  Vx  =  16  m/s) 


There  was  no  significant  variation  between  k-s  and  RSM  model  solutions,  except 
in  prediction  of  vortex  center  location  at  x/c  =  0.75  and  0.95.  The  k-s  model  better 
predicted  that  location  for  x/c  =  0.75,  and  RSM  model  better  predicted  it  for  x/c  =  0.95. 
This  was  likely  a  result  of  applying  averaging  and  steady  modeling  to  an  inherently 
unsteady  flow,  and  there  is  no  reason  to  believe  that  even  the  experimental  data  show  a 
“true”  vortex  core  position  since  breakdown  has  been  observed  to  shift  its  location 
(Morton,  Forsythe,  Mitchell,  Hajek,  2002;  Novak,  Sarpkaya,  1999;  Leibovich,  1978; 
Faler,  Leibovich,  1977). 
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Figure  3.42  -  Cp  Comparison  of  Steady  Turbulence  Models  at  x/c  =  0.55 
(Mesh  Revision  L,  a  =  18°,  Vx  =  16  m/s) 


Figure  3.43  -  Cp  Comparison  of  Steady  Turbulence  Models  at  x/c  =  0.75 
(Mesh  Revision  L,  a  =  18°,  Voo  =  16  m/s) 
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Figure  3.44  -  Cp  Comparison  of  Steady  Turbulence  Models  at  x/c  =  0.95 
(Mesh  Revision  L,  a  =  18°,  Ko  =  16  m/s) 


Figure  3.45  shows  nondimensional  velocity  contours,  defined  as  V/Vx,  contained 
in  the  plane  perpendicular  to  the  wing  surface  and  which  runs  through  the  center  of  the 
delta  wing’s  vortex.  Its  images  include  the  predicted  solution  from  the  steady  laminar,  S- 
A,  RNG  k-s,  and  RSM  flow  models.  In  the  figure,  predicted  vortex  breakdown  is 
indicated  where  the  vortex  diameter  expands  and  where  flow  may  reverse  and  stagnate, 
shown  by  regions  of  dark  blue,  where  flow  is  at  or  near  zero.  All  but  the  S-A  model 
predicted  stagnated  flow.  However,  experimentalists  have  observed  in  cases  where  Re  > 
3  x  l(f  (which  applies  to  this  case  where  Re  =  3.4  x  105)  that  in  the  region  of  vortex 
breakdown  and  immediately  following,  the  flow  does  not  necessarily  reverse,  stagnate  or 
bridge  between  laminar  and  turbulent  states  (Novak,  Sarpkaya,  1999:  825,  833).  Thus 
the  S-A  model  may  arguably  have  predicted  vortex  breakdown  but  not  as  pronounced  as 
by  the  other  models.  Additionally,  since  reversed  flow  results  in  bubble  type  breakdown, 
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Steady  RNG  k-c,  a  =  1 8°,  V^-  16  m/s  Steady  RSM.  a  =  1 8°,  Vin)  =  16  m/s 

Figure  3.45  -  Nondimensional  Velocity  Contours  through  Vortex  Center,  Predicted  by 
Steady  Turbulence  Models  (Mesh  Revision  L,  a  =  18°,  Vx  =  16  m/s) 

shown  in  Figures  1.3a  and  1.4,  the  vortex  core  diameter  naturally  becomes  large  -  as 

predicted  by  the  laminar,  k-s  and  RSM  models,  and  as  indicated  in  Figures  3.41-45  where 

these  models  predicted  a  larger  region  of  vortex  bursting  and  turbulence.  But  since 

experimental  findings  indicate  that  flow  at  this  Reynolds  number  typically  results  in 

spiral  type  breakdown  with  smaller  vortex  diameter,  the  S-A  turbulence  model  most 

accurately  predicted  a  solution  which  correlates  with  experimental  numerical  data  and 

with  experimental  observation  of  the  flow  physics,  at  least  with  respect  to  vortex  size 
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after  bursting.  This  conclusion  must  include  the  caveat  that  the  steady  S-A  model’s 
predicted  solution  was  only  marginally  acceptable,  since  Cp  magnitudes  were  low  by  as 
much  as  21%  from  the  experimental  value  and  since  predicted  vortex  diameter  at  x/c  = 
0.35  and  0.55  was  larger  than  indicated  experimentally  -  shown  by  lesser  Cp  curve  slope 
in  Figures  3.41  and  3.42.  Furthermore,  the  LSU  study,  which  used  the  same  wing  but 
attached  a  fuselage,  showed  distinct  and  relatively  large  regions  of  stagnated  flow  as  the 
vortex  approached  the  wing’s  trailing  edge  (Guillot,  1999:  41).  To  summarize,  the  S-A 
model  may  become  the  most  accurate  turbulence  model  after  appropriate  numerical  code 
modifications  are  incorporated  to  help  it  better  predict  the  physics  of  vortex  breakdown 
(Morton,  Forsythe,  Mitchell,  Hajek,  2002). 

Another  observation  from  the  S-A  model’s  predicted  solution  in  Figure  3.45  is 
that  the  vortex  rotational  speed  prior  to  breakdown  was  approximately  2,200  rad/sec  or 
350  Hz,  based  on  radial  velocity  and  distance  from  vortex  center.  This  result  is  less  than 
it  should  be  because  the  radial  distance  was  over-predicted  and  is  considerably  less  than 
an  experimental  observation  of  1,000  Hz  for  a  similar  configuration  and  of  860  Hz  for  the 
LSU  study  (estimated  graphically,  not  measured  specifically)  and  serves  as  further 
evidence  that  the  numerical  model  is  not  entirely  accurate  (Novak,  Sarpkaya,  1999:  825; 
Guillot,  1999:  41). 

Figure  3.46  shows  blue  particle  streamlines  through  the  vortex  and  bursting 
region  and  a  red  vortex  core  centerline  for  the  same  turbulence  models  shown  in  Figure 
3.45.  The  vortex  core  centerline  was  created  using  FIELD  VIEW’S  Vorticity  Alignment 
method,  and  vortex  breakdown  follows  shortly  after  the  first  large  break  in  that  line 
(Intelligent  Light,  2001 :  137).  The  streamlines  indicate  roughly  where  vortex  breakdown 
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Steady  RNG  k-e,  a  =  18°,  V^  -  16  m/s  Steady  RSM.  a  =  18°,  V^  -  16  m/s 

Figure  3.46  -  Blue  Streamlines  through  Vortex  and  Red  Vortex  Core,  Predicted  by 
Steady  Turbulence  Models  (Mesh  Revision  L,  a  =  18°,  Vx  =  16  m/s) 

occurs,  shown  by  radial  expansion,  and  they  show  reversed  flow  in  the  laminar,  k-s  and 

RSM  cases.  In  the  S-A  case,  top  right  in  Figure  3.46,  vortex  bursting  is  indicated  by  the 

first  large  break  in  the  vortex  core  centerline;  the  streamlines  relax  in  frequency  rather 

than  expand  significantly  in  core  diameter.  It  is  debatable  whether  the  S-A  model  even 

predicted  vortex  breakdown. 

From  Figures  3.45  and  3.46,  the  approximate  location  may  be  identified  where 
vortex  breakdown  occurs.  It  is  seen  in  Figure  3.45  where  there  begins  to  be  velocity 
gradients  inside  the  vortex  core  and  in  Figure  3.46  where  streamlines  show  relaxed 
rotational  frequency  or  vortex  diameter  expansion  and  where  a  large  break  occurs  in  the 
vortex  core  centerline  prediction.  With  each  turbulence  model,  the  solution  appears  to 
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predict  a  vortex  breakdown  location  between  x/c  =  0.35  and  0.45.  A  more  precise 
location  would  be  arguable,  since  it  is  a  subjective  detennination. 

These  results  agree  relatively  closely  with  three  experimental  studies.  Two  of 
those  studies  used  a  delta  wing  with  the  same  aspect  ratio  and  sweep  angle,  a  slightly 
greater  thickness,  and  Re  =  1.0  x  105,  and  those  studies  for  a  =  18°  each  determined 
vortex  breakdown  occurred  between  x/c  =  0.42  and  0.45  (Johari,  Olinger,  Fitzpatrick, 
1995:  806;  O’Neil  et  al,  1989).  The  third  study  (LSU)  used  the  same  half  delta  wing  as 
this  study  but  with  fuselage  attached;  in  that  configuration  with  a  =  15°,  one  may  argue 
the  equivalent  angle  of  attack  was  about  18°,  since  there  is  close  CL  agreement  seen  in 
Table  3.2  (using  results  from  the  steady  S-A  model  at  a  =  18°).  The  LSU  study 
determined  that  vortex  breakdown  occurred  between  x/c  =  0.3  and  0.4  (Guillot,  1999: 
38).  Traub’s  method  predicts  vortex  breakdown  occurs  at  x/c  =  0.23  for  a  =  18°  and  at 
x/c  =  0.43  for  a  =  15°,  based  on  empirical  data  for  delta  wings  with  larger  sweep  angles 
than  that  used  in  this  study  (1996;  see  Appendix  B  for  calculations).  An  earlier  study 
showed  breakdown  at  about  x/c  =  0.6  for  a  =  18°,  but  that  experiment  used  a  delta  wing 
with  a  thickness  of  0.25  cm  (about  20%  the  thickness  used  in  this  study),  sharper  leading 
edge  bevel  and  a  rough  upper  wing  surface  (Wentz,  Kohlman,  1971). 

Figure  3.47  shows  contour  slices  of  total  pressure  over  the  wing  surface  from  the 
steady  laminar  model’s  predicted  solution.  The  slices  are  perpendicular  to  the  upper 
wing  surface,  spaced  from  x/c  =  0.05  to  1.05  by  increments  0.10c.  In  this  figure,  vortex 
breakdown  occurred  where  the  inner  core  pressure  contours  first  expand,  which  was 
possibly  as  soon  as  between  x/c  =  0.25  and  0.35  but  more  likely  between  x/c  =  0.35  and 
0.45.  This  is  a  visually  qualitative  assessment  but  has  potential  to  feed  into  a  quantitative 
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Figure  3.47  -  Contours  of  Total  Pressure  over  Wing  Surface  for  Steady,  Laminar  Model 

(Mesh  Revision  L,  a  =  18°,  V«,  =  16  m/s) 


feedback  control  system,  should  the  system  need  to  identify  where  vortex  breakdown 
occurs.  This  steady  laminar  model  resulted  in  over-prediction  of  vortex  extent  but  it 
adequately  demonstrates  the  physics  of  vortex  bursting. 

Figure  3.48  shows  contour  slices  of  velocity  magnitude  nondimensionalized  by 
the  freestream  velocity  over  the  wing  surface  from  the  steady,  S-A  turbulence  model’s 
predicted  solution  at  the  same  locations  described  for  Figure  3.47.  Contours  were  limited 
between  V/Vx  =  1.0  and  1.6  for  easier  visibility  of  multiple  slices,  and  a  dark  blue  line 
indicates  the  vortex  core  centerline  by  vorticity  alignment.  This  contour  plot  was 
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Figure  3.48  -  Contours  of  Nondimensional  Velocity  over  Wing  Surface  for  Steady,  S-A 
Turbulence  Model  (Mesh  Revision  L,  a  =  18°,  Vx  =  16  m/s) 


included  for  comparison  with  equivalent  plots  of  data  generated  in  the  LSU  experiments, 
and  these  contours  compare  reasonably  well  with  those  experimental  results.  Forms  are 
basically  the  same  but  the  numerically  predicted  magnitudes  are  lower  by  roughly  10%  in 
the  areas  of  maximum  velocity.  Interestingly  the  LSU  data  indicated  no  regions  of 
stagnated  axial  flow  until  x/c  =  0.75  (Guillot,  1999:  41).  This  indicates  that  the  steady  S- 
A  model  is  perhaps  the  least  inaccurate  of  the  turbulence  models  evaluated  in  this  study, 
since  the  other  models  (laminar,  k-s  and  RSM)  predicted  stagnation  as  early  as  x/c  = 
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0.50;  however,  the  S-A  model  is  far  from  accurate,  since  it  did  not  predict  any  stagnation 
whatever  and  since  it  questionably  even  predicted  vortex  breakdown. 

As  one  additional  argument  against  the  steady  S-A  model’s  ability  to  correctly 
predict  vortex  breakdown,  Figure  3.49  compares  iso-surface  helicity  contours  for  the 
steady  S-A  and  laminar  models.  Helicity  is  the  dot  product  of  vorticity,  which  is  a 
measure  of  rotational  fluid  flow,  and  the  velocity  vector  (FLUENT,  2001:  27.4).  This 
figure  shows  that  the  laminar  model  predicts  breakdown  of  this  relatively  uniform 
rotational  flow  much  sooner  chordwise  than  it  does  with  the  S-A  model.  Interestingly, 
neither  model  predicted  breakdown  of  the  secondary  vortex  along  the  leading  edge, 
which  is  contrary  to  experimental  findings  that  the  secondary  vortex’s  breakdown 
preceded  that  of  the  primary  vortex  (Cummings,  Morton,  Siegel,  2003). 


Steady  S-A,  a  =  18°,  K,=  16  m/s  Steady  Laminar,  a  =  18°,  K-=  16  m/s 

Figure  3.49  -  Iso-Surface  Helicity  Contours  for  Steady  S-A  and  Laminar  Models 
(Mesh  Revision  L,  a  =  18°,  Vm  =  16  m/s) 


To  summarize  the  steady  solutions  from  the  various  non- turbulence  and 
turbulence  models,  the  inviscid,  laminar,  RNG  k-s,  and  RSM  models  over-predicted  the 
expanse  of  vortex  breakdown  and  its  ensuing  turbulence  and  under-predicted  vortex 
strength,  indicated  by  lower  velocity  and  higher  static  pressure  around  the  primary 
vortex.  The  steady  S-A  turbulence  model  also  under-predicted  vortex  strength,  though 
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generally  not  as  poorly  as  the  other  models;  except  at  x/c  =  0.35  (which  may  be  a  grid 
resolution  issue),  it  predicted  satisfactorily  the  region  and  extents  of  vortex  breakdown 
and  its  turbulent  wake.  Thus  none  of  the  steady  turbulence  models  predicted  a 
completely  satisfactory  solution. 

Since  “the  effect  of  grid  density  or  flow  solver  may  have  a  higher-order  effect  on 
the  solution  than  the  turbulence  model  itself’  (Krai,  1998:  485),  a  study  of  grid  resolution 
follows. 

Mesh  Adaptation 

As  part  of  a  mesh  independence  study  and  in  hopes  of  improving  the  numerically 
predicted  solutions,  Mesh  Revision  L  was  adapted  in  two  ways  using  the  Mesh 
Adaptation  features  in  FLUENT,  and  it  was  modified  also  in  Gridgen  to  form  Revision 
N.  The  first  adaptation  was  motivated  from  a  numerical  study  analyzing  delta  wing 
vortex  generation  with  an  unstructured  grid,  which  concluded  that  topological  features 
may  be  effectively  used  for  mesh  refinement  to  obtain  better  vortex  resolution  and  stated 
that  “grid  resolution  around  the  vortex  core  is  very  important  for  the  accurate  prediction 
of  the  vortex  breakdown”  (Murayama,  Nakahashi,  Sawada,  2001 :  1305,  131 1).  Also,  it  is 
necessary  to  use  “sufficient  grid  resolution  in  regions  of  high  flow  gradients  [such  as 
around  and  within  the  vortex  core]  to  obtain  accurate  numerical  solutions”  (Ekaterinaris, 
Schiff,  1990:  60).  In  accordance  with  these  recommendations,  a  FLUENT  Iso-Value 
Adaption  was  performed  on  the  steady  S-A  solution  from  Mesh  Revision  L  with  a  =  18° 
and  Voo  =  16  m/s.  This  adaptation  was  focused  on  iso-metric  contours  of  turbulent 
viscosity,  since  these  contours  captured  the  turbulent  field  within  the  vortex  core;  all  cells 
with  turbulent  viscosity  values  between  0.0019  and  0.0020  m /s  were  then  refined. 
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FLUENT  refines  three-dimensional  cells  by  dividing  each  respective  face  into  four  equal 
triangles  or  squares  then  propagating  it  volumetrically  (FLUENT,  2001:  23.2.2).  The 
result  of  this  refinement  is  shown  in  Figure  3.50,  where  the  image  on  top  shows  a  vertical 


Figure  3.50  -  Mesh  Adaptation  within  Vortex  Core 
(Revision  L,  Steady  S-A,  a  =  18°,  F®  =  16  m/s) 


RevL  (a=  18i de'g).  (S-A  steadi 

mesh  adapted  by  modified  Turbulent  viscosity  contour 
vertical  slice'  through  wing-  at  x/c=0.45 
transparent  nu_t  =  0.002 
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slice  through  the  turbulent  viscosity  three-dimensional  field  and  the  image  on  the  bottom 
shows  the  extent  of  the  adaptation  via  a  slice  through  the  vortex  core.  This  adaptation 
created  an  additional  327,000  tetrahedral  cells  for  a  56%  increase  in  total  number  of 
volumetric  cells. 

Another  adaptation,  independent  of  the  adaptation  described  above,  was 
performed  to  better  show  generation  of  the  secondary  vortex  along  the  wing’s  leading 
edge  and  incidentally  improved  portions  of  the  solution  prediction.  This  adaptation  used 
FLUENT’S  Region  Adaption  menu  option  and  created  a  cylindrical  region  for  refinement. 
The  cylinder’s  centerline  was  located  parallel  to,  1  cm  above  and  1  cm  inboard  of  the 
wing  leading  edge;  cylinder  radius  was  2.5  cm.  This  allowed  for  refined  cells  along  both 
lower  and  upper  surfaces  of  the  wing  around  the  leading  edge.  This  improved  resolution 
in  the  region  of  the  secondary  vortex  but  also  in  the  region  where  the  primary  vortex  is 
initialized  (or  where  the  shear  layers  separate),  resulting  in  a  generally  improved 
predicted  solution.  This  refinement  is  shown  in  Figure  3.51  as  a  tight  clustering  of  cells 
on  the  wing  surface  (bottom  image)  and  in  a  vertical  slice  around  the  wing  (top  image). 
An  additional  295,000  tetrahedral  cells  resulted  from  this  adaptation. 

Mesh  Revision  N  was  created  in  an  effort  to  refine  the  number  of  cells  around  the 
entire  wing  upper  surface  but  not  to  the  degree  generated  by  FLUENT.  A  three- 
dimensional  permeable  enclosure  was  created  with  roughly  a  2.5-cm  gap  between  it  and 
the  wing,  excluding  the  wing’s  bottom  surface;  it  resembles  a  box  for  a  slice  of  pie  or 
pizza,  as  shown  in  the  top  image  in  Figure  3.52.  The  bottom  image  shows  the  cell  face 
resolution  and  nodal  dimensions  on  the  upper  wing  surface,  which  dimensions  are  double 
those  of  Revision  L’s  upper  wing  surface.  These  changes  resulted  in  an  additional 
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Figure  3.5 1  -  Mesh  Adaptation  along  Leading  Edge 
(Revision  L,  Steady  S-A,  a=  18°,  F®  =  16  m/s) 
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Adapted  over  Upper  Surface,  Mesh  Revision  N  Adapted  over  Upper  Surface.  Mesh  Revision  N 

Figure  3.52  -  Mesh  Adaptation:  Revision  N  (Steady  S-A,  a  =  18°,  ¥«,  =  16  m/s) 


215,000  tetrahedral  cells  and  gave  results  comparable  to  those  from  the  leading-edge 
adaptation  but  at  less  computational  cost.  All  of  the  adaptations  lowered  y+  values  along 
the  wing  surface  but  still  not  sufficiently  to  meet  the  resolution  requirement  for  the  LES 
model. 
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Each  adaptation  was  used  to  generate  a  numerical  solution  using  S-A  turbulence 
model  and  steady  solver  with  a  =  18°  and  Voo  =  16  m/s.  Table  3.3  summarizes  computed 
lift  coefficients,  and  additional  comparisons  follow.  Solutions  from  the  adapted  meshes 
resulted  in  no  more  than  0.6%  variation  from  the  non-adapted  baseline’s  Cl,  thus 
showing  that  the  initial  numerical  mesh  allowed  for  adequate  prediction  of  this  integrated 
quantity. 

Table  3.3  -  Comparison  of  Cl  Resulting  from  Various  Mesh  Adaptations  and  Using 
Steady  S-A  Model  with  a  =  18°  and  14,  =  16  m/s 


Case 

CL 

Rev  F,  No  Adaptation 

0.7974 

Rev  F,  Turbulent  Viscosity  Adaptation 

0.8022 

Rev  F,  Feading  Edge  Region  Adaptation 

0.7978 

Rev  N,  Pie  Slice  Box  Adaptation 

0.7959 

The  right-side  image  in  Figure  3.53  shows  noticeable  improvement  in  contour 
smoothness  for  the  model  with  refined  mesh  within  the  vortex  core.  Smoother  contours 


for  this  adapted  mesh  are  also  visible  in  the  top  right  image  in  Figure  3.55.  In  addition  to 


finer  resolution,  this  mesh  adaptation  resulted  in  improved  Cp  slopes  and  magnitudes  in 


RevL  ( a  =  1 8  drs).  .S-A  steady 

Nondimcnsional  Velocity  Contour*  (u_mf«16  05  tn/s) 
mesh  Adapted  by  modified  turbulent  viscosity  contours 
longitudinal^tc^-titong  vurt**-c«£e  centerline 


No  Adaptation,  Mesh  Revision  L  Adapted  within  Vortex,  Mesh  Revision  L 

Figure  3.53  -  Nondimensional  Velocity  within  Vortex  Core  for  Revision  F  and 
Adaptation  by  Turbulent  Viscosity  Iso-Contours  (Steady  S-A,  a  =  18°,  F®  =  16  m/s) 
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Figures  3.56-59;  its  prediction  (blue  line  with  no  markers)  exceed  the  non-adapted 
solution  (green  line)  by  4-6%.  However,  these  improvements  were  small  and  were  not 
deemed  worthwhile  for  the  tradeoff  in  computational  speed. 

Figure  3.54  compares  plots  of  no-slip  surface-restricted  flow  (sometimes  termed 
oil  surface  flow)  predicted  by  the  steady  S-A  models  with  mesh  adaptation  along  the 
leading  edge  (top  right),  over  the  upper  wing  surface  (bottom),  and  with  the  non-adapted 
mesh  (top  left).  In  addition  to  providing  greater  resolution  of  the  surface  flow  lines,  both 
adapted  mesh  solutions  better  indicate  presence  of  a  secondary  vortex,  along  the  leading 
edge.  Surface-restricted  flow  plots  showed  little  variation  among  solutions  from  different 


No  Adaptation,  Mesh  Revision  L 


Adapted  along  Leading  Edge.  Mesh  Revision  L 


Adapted  over  Upper  Surface,  Mesh  Revision  N 

Figure  3.54  -  Wing  Surface-Restricted  Flow  for  Mesh  Revision  L,  Adaptation  along 
Leading  Edge,  and  Revision  N  (Steady  S-A,  a  =  18°,  Vx=  16  m/s) 
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turbulence  models,  as  well  as  between  steady  and  unsteady  solutions;  thus  no  additional 
oil  flow  plots  are  included  in  this  document. 

Improved  secondary  vortex  resolution  becomes  more  apparent  after  observing 
total  pressure  contours  in  Figure  3.55.  In  that  figure’s  images  from  the  adapted  mesh 
solutions  (bottom  left  and  right),  for  adaptation  along  the  leading  edge  and  over  the  upper 
surface,  the  secondary  vortex  along  the  leading  edge  has  greater  strength  and  is  better 
distinguished  from  the  primary  vortex.  It  appears  further  that  the  leading  edge  mesh 


No  Adaptation,  Mesh  Revision  L 


Adapted  within  Vortex,  Mesh  Revision  L 


WevL  l«  !B  J^s!  B  A  slrndv - 

Total  Pressure  Contours  (Pa) 

Fnesh  adapted  by  modified  turbulent  viscosity  contours 
Rlices  from  Ob  to  I  Ob  by  10  increments 


Total  Pressure  Contours  (Pa) 

RevL  (a-  18  deg).  S-A  steady 

llieM  from  x/c-  05  to  105  by  10  increments 


rolal  Wrssure  rWours  U’aJ - 

RevL  (a-18  deg).  S-A  steady 
mesh  adapted  along  2ndory’  vortex 
slice*  from  05  to  I  05  by  10  increments 


Kvv\  1H  drgl.  •'  A  si  rad  v - 

total  Pressure  Contours  (Pn) 

mesh  refined  over  upper  wing  surface 

slices  from  05  to  1  05  by  10  increments 


Adapted  along  Leading  Edge,  Mesh  Revision  L 


Adapted  over  Upper  Wing  Surface,  Mesh  Revision  N 


Figure  3.55  -  Contours  of  Total  Pressure  over  Wing  Surface  for  Mesh  Revision  L, 
Adaptations,  and  Revision  N  (Steady  S-A,  a  =  18°,  V«,  =  16  m/s) 
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refinement  gave  better  resolution  than  did  Mesh  Revision  N,  indicated  also  in  Figure 
3.56,  where  a  small  secondary  Cp  peak  is  evident  around  y/s  =  0.34;  but  again  it  is  an 
issue  of  computational  load,  where  Revision  N  iterated  and  converged  more  quickly. 

Apparently,  mesh  adaptation  using  turbulent  viscosity  contours  failed  to  refine  the 
vortex  core  in  its  entirety,  as  the  other  adaptations  appear  to  have  better  refined  the  vortex 
core  centerline,  indicated  by  a  thick  blue  line  in  each  image  in  Figure  3.55.  Nonetheless, 
each  adapted  case  resulted  in  improved  centerline  resolution,  which  may  assist  in 
identification  of  vortex  breakdown  location. 

Figure  3.56  shows  quite  clearly  that  mesh  adaptation  improved  the  solution  at  and 
likely  around  x/c  =  0.35;  this  shows  also  that  the  solution  from  Revision  L  was  not  mesh 
independent  at  that  location.  Comparing  against  the  non-adapted  results,  the  numerical 
model  with  leading  edge  adaptation  improved  the  peak  Cp  magnitude  by  13%,  and  the 
two  other  adapted  meshes  each  resulted  in  5.7%  magnitude  increase.  They  also  slightly 
improved  the  CP  curves’  slopes  or  vortex  extents.  However,  Figures  3.57-59  indicate  that 
Revision  L’s  solution  was  independent  of  mesh  at  most  other  locations  over  the  wing 
surface.  Therefore,  since  the  regions  corresponding  to  x/c  >  0.35  have  greater 
significance  for  vortex  breakdown  control  issues,  the  non-adapted  Mesh  Revision  L  was 
considered  sufficient  for  further  evaluation.  As  an  aside,  the  model  used  for  follow-on 
investigations  would  do  well  to  include  improved  grid  resolution  in  the  region 
corresponding  to  the  forward  half  of  the  wing. 
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Figure  3.56  -  Cp  Comparison  of  Mesh  Revision  N  and  Adaptations  of  Mesh  Revision  L 
at  x/c  =  0.35  (Steady  S-A,  a  =  18°,  Vm  =  16  m/s) 


Adapted  Grids  (Steady  S-A,  a=18  deg,  u_inf=16  m/s) 
x/c  =  0.55 


Figure  3.57  -  Cp  Comparison  of  Mesh  Revision  N  and  Adaptations  of  Mesh  Revision  L 
at  x/c  =  0.55  (Steady  S-A,  a  =  18°,  Vm  =  16  m/s) 
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Figure  3.58  -  Cp  Comparison  of  Mesh  Revision  N  and  Adaptations  of  Mesh  Revision  L 
at  x/c  =  0.75  (Steady  S-A,  a  =  18°,  F®  =  16  m/s) 


Figure  3.59  -  Cp  Comparison  of  Mesh  Revision  N  and  Adaptations  of  Mesh  Revision  L 
at  x/c  =  0.95  (Steady  S-A,  a  =  18°,  F®  =  16  m/s) 
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In  retrospect,  it  was  a  poor  decision  to  not  proceed  at  least  with  Mesh  Revision  N, 
because  it  was  later  learned  that  the  secondary  vortex  breakdown  appears  to  affect  and 
perhaps  feed  into  the  primary  vortex  breakdown  (Cummings,  Morton,  Siegel,  2003). 
With  that  revelation,  it  would  be  prudent  to  use  a  model  which  predicts  and  adequately 
resolves  the  secondary  vortex  flow  properties,  as  well  as  those  of  the  primary  vortex. 
Whereas  in  this  study,  improved  resolution  of  the  secondary  vortex  was  deemed  an 
academically  interesting  but  non-essential  exercise,  and  its  resultant  improvement  in  the 
overall  solution  was  weighed  unfavorably  against  the  required  additional  computation 
time  and  resources. 

Effect  of  Turbulence  Model  with  Time  Accurate  Flow 

This  section  discusses  some  pros  and  cons  of  using  an  unsteady  and  hence  more 
accurate  solver  for  this  highly  unsteady  problem  and  comparison  of  the  S-A,  RNG  k-s, 
and  RSM  turbulence  models  using  Mesh  Revision  L  with  a  =  18°  and  Vx  =  16  m/s. 
Inviscid  and  laminar  models  were  not  evaluated  in  the  unsteady  domain  due  to  the 
turbulent  nature  of  the  flow  within  and  aft  of  vortex  breakdown.  While  the  LES  model 
(unsteady  by  default)  was  used  preliminarily,  its  results  and  poor  correlation  with 
experimental  data  confirmed  the  requirement  for  a  fine  numerical  mesh  (with  y+  ~  1 
along  all  surfaces  of  interest);  due  to  limitations  on  time  and  computing  resources  in  this 
study,  the  fine  computational  grid  was  not  completed.  Adapted  meshes  were  not 
evaluated  in  time  accurate  flow,  though  that  is  a  recommended  improvement  for  future 
work. 

Using  an  unsteady  solver  for  this  problem  had  certain  advantages  over  using  a 
steady  solver.  As  expected  it  predicted  flow  physics  which  were  considerably  more 
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consistent  with  theoretical  assertions  and  experimental  findings.  Once  flow  was  fully 
developed  (after  roughly  2  sec  in  dimensional  time),  Figure  3.60  shows  with  streamlines 
through  the  vortex  core,  using  the  RSM  turbulence  model,  that  there  was  notable 
fluctuation  over  time  in  the  flow  structure  within  and  following  vortex  breakdown  -  in 
radial  extent  of  the  core  and  in  local  flow  direction  within  the  core.  A  steady  solver  has 


Figure  3.60  -  Streamlines  through  Vortex  from  Unsteady  RSM  Turbulence  Model 

(Revision  L,  a  =  1 8°,  Vx  =  1 6  m/s) 


not  the  capacity  to  predict  temporal  fluctuations  in  flow  structure.  However,  Figure  3.61 
shows  with  streamlines  through  the  vortex  core,  using  the  k-s  turbulence  model,  that 
there  was  no  fluctuation  over  time  in  the  flow  structure  and  that  the  unsteady  solution  is 
essentially  equivalent  to  the  time  averaged  solution  (top  left  image  in  the  figure).  The  S- 
A  model  yielded  also  a  time  independent  vortex  structure  and  one  which  failed  to 
adequately  predict  vortex  breakdown,  just  as  it  failed  to  do  in  steady  state.  Its 
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t  =  1 .68  sec 


t  =  1 .76  sec 


Figure  3.61  -  Streamlines  through  Vortex  from  Unsteady  RNG  k-s  Turbulence  Model 

(Revision  L,  a  =  18°,  F®  =  16  m/s) 

nondimensional  velocity  contours  through  the  vortex  core  after  flow  was  fully  developed, 
not  included  here  for  brevity,  are  virtually  identical  to  those  in  the  top  right  image  in 
Figure  3.45.  These  unsteady  turbulence  models’  (S-A  and  k-s)  inaccurate  predictions 
disqualify  them  as  candidates  for  a  suitable  unsteady  solver  for  this  problem. 

Again  from  the  unsteady  RSM  model,  Figure  3.62  compares  Cp  values  at  eight 
different  time  steps  between  t  =  1.52  and  2.08  sec  and  shows  this  model  failed  to 
accurately  predict  Cp  values  comparable  to  the  experimental  data.  At  x/c  =  0.35  and  0.55, 
these  data  indicate  that  the  flow  was  fully  developed  since  Cp  magnitude  and  peak 
locations  did  not  change  with  time.  At  x/c  =  0.75,  it  appears  there  was  fluctuation  in 
vortex  strength,  though  not  to  a  significant  degree.  At  x/c  =  0.95,  which  has  an  exploded 
view  in  Figure  3.63,  there  is  cyclic  variation  in  vortex  strength  and  in  centerline  location, 
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Figure  3.62  -  Cp  Plots  from  Unsteady  RSM  Turbulence  Model  for  t  =  1.52-2.08  sec 

(Revision  L,  a  =  18°,  Vx  =  16  m/s) 


Steady  and  Unsteady  RSM  Turb  Model  with  Rev  L 
t  =  1 .52  -  2.08  sec  by  0.08  increments 
x/c  =  0.95 


Figure  3.63  -  Cp  Plots  from  Steady  and  Unsteady  RSM  Turbulence  Models  for  t  =  1.52- 
2.08  sec  at  x/c  =  0.95  (Revision  L,  a  =  18°,  Fa,  =  16  m/s) 
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where  the  vortex  center  wandered  between  y/s  =  0.53  and  0.70,  which  corresponds  to  a 
spanwise  range  of  3.4  cm.  This  is  consistent  with  experimental  observations  that  the 
vortex  core  meanders  several  core  radii  in  any  given  direction.  However,  the  unsteady 
RSM  model  did  not  predict  this  core  meandering  farther  upstream  over  the  wing. 

Inspection  of  Figure  3.63  also  reveals  that  the  vortex  core  meanders  in  a  cyclic 
pattern.  Following  Cp  curves  by  time  history  on  the  figure,  they  translate  in  an  ovular, 
clockwise  direction  with  a  period  of  about  0.48  sec  or  frequency  of  roughly  2  Hz.  This  is 
consistent  with  experimental  observations  that  the  vortical  structure  rotates  at  1-2  Hz. 
This  vortical  rotation  is  further  evidenced  by  cyclic  variation  in  Cp  magnitude  in  Figure 
3.63;  when  the  vortex  core  rotates  away  from  the  wing  surface,  static  pressure  on  the 
surface  increases  and  results  in  a  lower  Cp  value.  Therefore  these  data  show  the  vortex 
core  rotates  outboard  and  upward  over  time  or  in  a  counter-clockwise  rotation,  which  is 
consistent  with  experimental  observation  that  the  vortex  structure  rotates  in  the  same 
direction  with  the  vortex  proper  -  both  are  counter-clockwise  for  the  port-side  wing  when 
looking  downstream.  Since  data  were  written  every  0.08  seconds  (200  time  steps  or 
4,000  computational  iterations),  it  was  not  possible  to  show  graphically  the  vortex 
rotational  frequency  of  several  hundred  cycles  per  second.  None  of  this  validation  with 
experimental  findings  was  possible  using  a  steady  solver.  In  summary,  the  unsteady 
solver  is  advantageous  over  steady  primarily  in  predicting  several  of  the  dynamic 
characteristics  of  the  flow.  It  appears  also  that  the  steady  RSM  model  predicted  the 
correct  time  averaged  vortex  core  centerline  position. 

Predicting  dynamic  characteristics  of  the  flow,  however,  may  or  may  not  be 
important  for  a  given  application.  The  main  disadvantages  of  using  an  unsteady  solver 
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are  the  time  required  to  achieve  fully  developed  flow  and  the  arguably  insignificant 
improvement  in  results.  An  unsteady  solution  is  obtained  in  weeks  rather  than  the  hours 
required  for  a  steady  solution.  There  is  no  question  whether  time  accurate  solutions, 
using  the  RSM  turbulence  model,  more  correctly  predict  the  unsteady  phenomena  of 
vortical  flow  and  vortex  breakdown,  but  is  prediction  of  these  occurrences  critical  to  the 
control  system  which  may  later  be  developed  with  assistance  from  a  numerical  flow 
model  and  solver?  This  question  of  course  needs  to  be  addressed  by  those  who  will 
develop  the  control  system. 

None  of  the  evaluated  unsteady  turbulence  models  were  able  to  predict  any 
variation  in  vortex  breakdown  location,  which  was  observed  experimentally.  Figures 
3.64-65  show  that  flow  was  fully  developed  not  later  than  about  t  =  1.6  sec  for  the 


unsteady  S-A  and  k-s  models,  since  the  predicted  CP  values  no  longer  appreciably 


Steady  and  Unsteady  S-ATurb  Model  with  Rev  L 
t=  1 .44  -  2.00  sec  by  0.08  increments 
left  to  right  -  x/c  =  0.35,  0.55,  0.75,0.95 


Figure  3.64  -  Cp  Plots  from  Unsteady  S-A  Turbulence  Model  for  t 

(Revision  L,  a  =  18°,  V«,=  16  m/s) 


1.44-2.00  sec 
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Steady  and  Unsteady  k-epsilon  Turb  Model  with  Rev  L  ^ 

t  =  1 .44  -  2.00  sec  by  0.08  increments 
lento  right -X(C=  0.35,  0.55,  0.75,0.95 


v _ / 

Figure  3.65  -  Cp  Plots  from  Unsteady  RNG  k-s  Turbulence  Model  for  t  =  1.44-2.00  sec 

(Revision  L,  a  =  18°,  V«>=  16  m/s) 


changed  with  time.  Different  from  the  RSM  predicted  solution,  these  models’  solutions 
showed  no  fluctuation  in  vortex  core  location  over  time.  This  inability  to  predict  the 
unsteady  flow’s  physics  essentially  negates  any  advantage  obtained  by  using  an  unsteady 
solver  for  these  two  turbulence  models.  In  fact,  unsteady  solutions  with  the  S-A  or  k-s 
turbulence  model  give  only  the  great  disadvantage  of  requiring  significantly  more 
computation  time  and  resources  than  the  corresponding  steady  solvers.  The  RSM 
turbulence  model  resulted  in  a  solution  with  cyclic  rotation  of  the  vortex  core  centerline 
but  aft  of  vortex  breakdown.  This  result,  seen  for  the  RSM  case  in  Figures  3.62,  3.63  and 
3.66,  is  good  and  bad  -  good  because  it  may  be  easier  to  control  breakdown  if  it  is 
predicted  to  remain  in  a  steady  location,  and  bad  because  it  failed  to  accurately  model  the 
physical  happenings  within  the  bursting  vortical  flow  field. 
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The  steady  RSM  solution  predicted  location  of  vortex  breakdown  just  as  well  as 
the  time  accurate  solver,  where  both  solvers  showed  its  occurrence  at  x/c  =  0.45.  Figure 
3.63’s  dark  blue  line  with  diamond  markers  shows  the  steady  RSM  prediction,  which 
appears  to  be  an  excellent  averaged  value  of  the  unsteady  predicted  Cp  values,  though  it 
does  not  show  local  extrema.  Figures  3.64-65,  which  also  indicate  the  steady  Cp  curves 
via  dark  blue  lines  with  diamond  markers,  show  that  the  steady  S-A  and  k-s  models 
adequately  averaged  the  unsteady  flow  properties. 

Total  pressure  contours  in  Figure  3.66  show  that  the  RSM  steady  solution  (top 
left)  closely  matches  any  of  the  unsteady  snapshots  of  the  solution,  where  the  small 
variations  occur  along  the  secondary  vortex  and  aft  of  x/c  =  0.65  for  the  primary  vortex. 
Nondimensional  velocity  contours  in  Figure  3.67  are  quite  similar  for  steady  and 
unsteady  RSM  predictions  everywhere  except  within  the  vortex  core  between  x/c  =  0.70 
and  1.20,  which  may  be  a  region  of  little  significance  for  vortex  breakdown  control. 
Figure  3.68  shows  nondimensional  velocity  contours  predicted  by  k-s  steady  and 
unsteady  solvers  and  indicates  that  both  predict  the  same  basic  vortex  and  vortex 
breakdown  structures.  Figures  3.67  and  3.68  further  indicate  that  neither  the  S-A  nor  k-s 
models  predicted  any  fluctuation  in  vortex  core  location  over  time  and  that  their  steady- 
state  solutions  were  essentially  equivalent.  It  is  worth  repeating  that  none  of  the 
turbulence  models  evaluated  in  this  study  adequately  predicted  wing  surface  pressure 
values  and  vortex  strength  and  size,  when  compared  to  the  UC  experimental  data. 
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Figure  3.66  -  Total  Pressure  Contours  from  Steady  and  Unsteady  RSM  Turbulence 
Models  (Revision  L,  a  =  18°,  Too  =  16  m/s) 
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Figure  3.67  -  Nondimensional  Velocity  Contours  through  Vortex  Core  from  Steady  and 
Unsteady  RSM  Turbulence  Models  (Revision  L,  a  =  18°,  Uo  =  16  m/s) 
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Figure  3.68  -  Nondimensional  Velocity  Contours  through  Vortex  Core  from  Steady  and 
Unsteady  RNG  k-s  Turbulence  Models  (Revision  L,  a  =  18°,  Vx  =  16  m/s) 


Figure  3.69  was  included  primarily  to  show  that  computational  solutions  may  be 
closely  compared  against  experimentally  determined  velocity  contours  within  the  vortical 
flow  field  (Guillot,  1999:  41).  This  figure  also  shows  that  the  vortex  was  steady  over 
time  at  x/c  =  0.35  but  that  its  higher  pressure  core  varied  considerably  over  time  at  x/c  = 
0.95. 
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Figure  3.69  -  Nondimensional  Velocity  Contours  at  x/c  =  0.35  and  0.95  from  Unsteady 
RSM  Turbulence  Model  (Revision  L,  a  =  18°,  F*  =  16  m/s) 
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Comparing  perfonnance  of  the  three  unsteady  turbulence  models  -  S-A,  RNG  k-s 
and  RSM,  Table  3.4  shows  the  unsteady  S-A  model  predicted  a  Cl  value  closest  to  the 
experimental  value,  while  the  k-s  and  RSM  predictions  were  within  6-8%  of  that  value. 
The  S-A  and  RSM  unsteady  solvers  predicted  Cl  values  which  improved  upon  the  steady 
prediction  by  1.0  and  1.6%,  respectively,  while  the  k-s  prediction  varied  only  by  0.1% 
between  steady  and  unsteady  solvers. 

Table  3.4  -  Comparison  of  Cl  for  Steady  and  Unsteady  (at  t  =  2.0  sec)  Turbulence 
Models  (Revision  L,  a  =  18°,  Vx  =  16  m/s) 


Turbulence 

Model 

CL 

steady  unsteady 

S-A 

0.797 

0.805 

RNG  k-s 

0.738 

0.737  i 

RSM 

0.740 

0.752 

Experiment 

0.8 

- 

Figure  3.70  shows  predicted  Cp  values  at  the  four  chordwise  locations  at  t  =  2.00 
sec,  where  flow  was  fully  developed  in  each  case.  These  results  are  quite  similar  to 
steady  results  seen  in  Figures  3.41-44,  effectively  casting  doubt  upon  the  unsteady 
solvers’  utility  for  this  problem.  At  each  chordwise  location  plotted  in  Figure  3.70,  the 
unsteady  S-A  turbulence  model  most  closely  predicted  peak  experiment  Cp  values  (or 
vortex  strength),  Cp  curve  slopes  (or  extent  of  the  generated  vortex),  and  peak  Cp  location 
(or  location  of  vortex  core  centerline).  Since  the  unsteady  S-A  model  failed  to  predict 
vortex  breakdown  or  any  unsteadiness  within  the  vortex  structure,  however,  it  is  an 
unacceptable  numerical  model.  Further,  based  on  results  in  Figure  3.70,  neither  the 
unsteady  k-s  nor  RSM  models  predicted  acceptable  pressure  values  over  the  wing 
surface,  excluding  them  also  as  viably  acceptable  numerical  models. 
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Figure  3.70  -  Cp  Plots  from  Unsteady  Turbulence  Models  at  t  =  2.00  sec 
(Revision  L,  a  =1 8°,  Va,  =  1 6  m/s) 


In  summary,  as  was  the  case  with  the  steady  turbulence  models,  none  of  the  time 
accurate  models  adequately  predicted  vortex  strength  or  dimensional  extents,  when 
compared  with  testing  data  from  UC.  Each  of  them  under-predicted  strength  of  the 
vortex  and  over-predicted  its  size.  While  the  unsteady  S-A  model  again  came  closest  in 
predicting  Cp  and  Cl  values,  it  failed  to  unquestionably  predict  occurrence  of  vortex 
breakdown  and  certainly  failed  to  correctly  predict  pressure  and  velocity  gradients  within 
the  vortex.  The  unsteady  RSM  turbulence  model  predicted  cyclic  meandering  of  the 
vortex  core  position,  though  not  as  far  upstream  as  was  noted  experimentally,  and 
predicted  no  variation  in  vortex  breakdown  location.  Unsteady  S-A  and  RNG  k-s  models 
categorically  failed  to  predict  any  of  the  unsteady  characteristics  of  vortical  flow  and 
vortex  breakdown  over  a  delta  wing. 
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There  are  no  clear  answers  to  whether  a  steady  model  may  suffice  for  this  control 
problem;  it  will  depend  largely  on  user  requirements,  which  as  yet  are  unspecified.  One 
possible  course  of  action  would  be  to  use  an  unsteady  solver  with  a  given  turbulence 
model  to  validate  accurate  correlation  with  the  flow  physics  then  use  its  corresponding 
steady  (quicker)  solver  to  predict  vortex  breakdown  position  given  various  initial 
conditions. 
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IV.  Numerical  Simulation  with  Flow  Control 


This  chapter  presents  results  for  along-core  blowing  to  control  vortex  breakdown 
over  the  half  delta  wing  numerical  flow  model.  In  the  following  sections,  the  numerical 
mesh  generation  and  adaptation  as  well  as  the  selected  flow  solver  parameters  will  be 
discussed  before  presenting  results. 

Mesh  Generation  and  Adaptation 

While  some  along-core  blowing  trial  runs  were  perfonned  on  earlier  grid  models, 
the  final  configuration  utilized  Mesh  Revision  L,  discussed  in  Chapter  III.  However,  it  is 
worthwhile  to  include  discussion  of  a  modification  to  the  initial  Revision  L,  a  mesh 
adaptation  of  Revision  L’s  blowing  axis,  and  creation  of  Revision  M. 

All  of  the  initial  meshes  included  the  wing’s  three  blowing  ports,  though  they 
were  set  as  a  wall  boundary  condition;  thus  they  served  both  as  placeholders  for  later 
blowing  and  as  markers  to  verify  that  predicted  location  of  the  vortex  core  centerline 
aligned  with  the  experimentally  identified  centerline.  From  these  early  solutions,  it 
appeared  the  vortex  core  was  not  properly  aligned,  so  when  one  of  the  ports  was  activated 
to  provide  blowing,  it  affected  neither  the  vortex  nor  its  breakdown.  This  result  agreed 
with  an  observation  in  the  experimental  LSU  study  that  when  blowing  is  angled  to  miss 
the  vortex  core,  the  jet’s  mass  flow  joins  the  larger,  swirling  vortex  flow  and  does  not 
affect  breakdown  (Guillot,  1999). 

Before  moving  port  locations  to  better  align  with  the  vortex  core  to  generate 
qualitative  results  -  as  experimental  configuration  would  have  been  violated,  negating 
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quantitative  assessment  -  a  close  comparison  was  made  between  wing  drawing,  which 
was  not  the  fabrication  drawing,  and  a  digital  image  of  the  actual  wing.  Since  they  did 
not  agree,  the  ports  in  grid  Revision  L  were  moved  to  locations  determined  from  the 
image;  consequently,  these  locations  proved  to  be  aligned  with  the  vortex  core  centerline. 
The  ports  were  numbered  such  that  port  1  was  closest  to  the  wing  apex,  located  at  x/c  = 
0.30  and  y/s  =  0.209  (where  5  is  the  wing  half  span  or  19.8  cm);  port  2  was  located  at  x/c 
=  0.60  and  y/s  =  0.409;  and  port  3  was  closest  to  the  trailing  edge,  located  at  x/c  =  0.80 
and  y/s  =  0.543.  The  spanwise  locations  would  not  be  valid  for  different  wing  angles  of 
attack,  though  jet  angles  could  be  appropriately  modified  to  compensate  for  a  different 
vortex  location. 

With  this  corrected  configuration,  blowing  still  failed  to  affect  vortex  breakdown. 
It  was  not  discovered  until  later  that  this  failure  occurred  in  part  because  the  degree  of 
vortex  breakdown  was  not  significant  enough  to  merit  affectation.  That  is,  the  S-A 
turbulence  model  failed  to  sufficiently  predict  vortex  bursting,  or  pressure  and  velocity 
gradients  within  the  core,  for  the  blowing  to  take  effect  any  more  than  it  would  affect  a 
vortex  with  no  breakdown.  Nonetheless  at  that  time,  in  hopes  of  better  resolving  jet 
blowing  into  the  vortex  core.  Mesh  Revision  L  was  adapted  in  FLUENT  by  refining  cells 
in  the  region  defined  by  a  cylinder  0.5  cm  in  radius  and  6.0  cm  in  length,  where  the 
leading  end  or  cylinder  top  was  centered  on  the  port  and  aligned  normal  to  the  blowing 
vector,  as  shown  in  Figure  4.1.  Adaptation  was  performed  twice  on  this  region  and 
created  an  additional  1 1,000  nodes  and  56,000  tetrahedral  cells. 

After  the  adapted  mesh  did  not  improve  results,  still  before  realization  of  error  in 
turbulence  model  selection,  and  after  noting  minute  mass  flow  loss  at  the  port,  Revision 
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M  modified  the  numerical  mesh  to  allow  for  blowing  normal  to  the  port  surface.  Figure 
4.2  shows  a  shaded  visualization  of  the  grid  at  and  around  a  blowing  port;  this  revision 
modeled  more  precisely  the  test  article’s  blowing  ports.  Building  this  numerical  model 
required  great  attention  to  ensure  the  port  was  set  perfectly  normal  to  the  three- 


Figure  4.1  -  Mesh  Adaptation  to  Revision  L  in  Jet  Blowing  Region 


Figure  4.2  -  Mesh  Generation,  Revision  M:  Views  of  Recessed,  Angled  Blowing  Port 
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dimensional  blowing  vector.  Since  the  final  model  from  this  study  may  ultimately  be 
used  to  show  vortex  breakdown  sensitivity  to  different  blowing  angles,  this  revision  was 
abandoned  due  to  the  need  to  carefully  rebuild  a  grid  for  each  variation  of  blowing  angle. 
On  the  other  hand  for  Revision  L,  the  user  need  only  alter  the  directional  components  for 
each  blowing  port  boundary  condition.  The  final  answer  was  that  mesh  Revision  L  was 
suitable  for  the  blowing  cases.  However,  due  to  apparent  dissipation  of  the  jet  momenta, 
as  shown  later  in  this  chapter,  it  would  be  prudent  to  investigate  further  refinement  of  the 
numerical  mesh’s  port  blowing  regions. 

Solver  Parameters 

For  blowing  cases,  the  same  final  model  configuration  was  used  except  the  active 
port’s  boundary  condition  was  changed  from  impermeable  wall  to  either  velocity  or 
pressure  inlet  in  the  FLUENT  define^boundary  conditions  menu.  Specifying  a  velocity 
inlet  was  preferred  because  it  required  no  translation  between  pressure  and  velocity 
values;  the  user  need  only  input  the  desired  velocity,  combined  with  directional 
components.  While  FLUENT  warns,  “This  [velocity  inlet]  boundary  condition  is 
intended  for  incompressible  flows,  and  its  use  in  compressible  flows  will  lead  to  a 
nonphysical  result  because  it  allows  stagnation  conditions  to  float  to  any  level” 
(FLUENT,  2001:  6.4),  the  user  may  circumvent  the  issue  of  uncontrolled  stagnation 
conditions  by  specifying  an  outflow  gauge  pressure  (99.34  kPa  for  Port  1)  which 
corresponded  to  the  experimentally  measured  pressure  near  the  blowing  port.  If  those 
data  were  unavailable,  the  user  would  likely  need  to  use  the  pressure  inlet  boundary 
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specification.  Additionally,  since  most  of  the  flow  within  the  numerical  wind  tunnel  was 
indeed  uncompressed,  using  the  incompressible-preferred  velocity  inlet  at  the  port  was 
not  a  poor  assumption. 

Other  required  inputs  at  the  blowing  port  boundary  included  temperature,  set  at 
298  K,  either  or  both  of  turbulent  viscosity  ratio  (TVR)  and  turbulence  intensity, 
depending  on  turbulence  model  used,  and  Cartesian  components  for  the  blowing  vector. 
TVR  was  set  at  1  and  turbulence  intensity  was  0.05%,  both  indicating  laminar  flow  at  the 
blowing  port.  Components  for  the  blowing  vector  were  calculated  using  the 
transformation  matrix  found  in  Appendix  D.  For  a  =  18°,  jet  elevation  angle  =  35°  and 
jet  azimuthal  angle  =  155°,  they  were  detennined  to  be  x-component  =  0.883313,  y- 
component  =  0.346189,  and  z-component  =  -0.3 16088.  The  resultant  blowing  vector  was 
angled  precisely  into  the  center  of  the  vortex  core.  All  steady  and  unsteady  blowing 
cases  were  partitioned  to  compute  in  parallel  on  between  12  and  18  processors. 

A  sub-study  was  performed  to  assess  any  adverse  ramifications  of  using  a  velocity 
inlet  instead  of  the  FLUENT-recommended  pressure  inlet  for  compressible  flow.  For 
each  case  the  steady  solver  with  S-A  turbulence  model  was  used,  with  a  =  18°  and  Coo  = 
16.05  m/s.  Comparisons  were  made  for  blowing  velocities  of  200  and  250  m/s.  For  the 
pressure  inlet  cases,  initial  pressure  was  set  at  99.56  kPa  and  total  gauge  pressures  were 
set  at  126.8  and  146.1  kPa,  corresponding  to  the  respective  velocities.  Comparison  of  Cl 
in  Table  4.1  shows  that  difference  between  using  pressure  or  velocity  inlet  is  negligible. 
Further,  comparison  of  Cp  in  Figure  4.3  shows  no  difference  between  the  two  inlet 
conditions  for  the  case  with  blowing  at  250  m/s;  the  case  with  blowing  at  200  m/s  had 
comparable  results. 
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Table  4.1  -  Comparison  of  Cl  from  Varied  Inlet  Specification  and  Port  Blowing  Velocity 


Case 

Viet  (m/s) 

CL 

Pressure  Inlet 
Velocity  Inlet 

200.0125 

200 

0.8054 

0.8091 

Pressure  Inlet 
Velocity  Inlet 

250.0017 

250 

0.8048 

0.8079 

Figure  4.3  -  Cp  Comparison  from  Velocity  and  Pressure  Inlet  at  Blowing  Port  1 
(Revision  L,  Steady  S-A,  a  =  18°,  Voo  =  16  m/s) 


Because  the  S-A  turbulence  model  did  not  adequately  predict  vortex  breakdown 
and  because  the  two-equation  RNG  k-s  model  was  lower  in  fidelity  by  definition,  the 
seven-equation,  time  accurate  RSM  turbulence  model  was  used  to  evaluate  effects  of 
along-core  blowing  into  the  vortex.  Since  this  turbulence  model  previously  demonstrated 
its  inability  to  accurately  predict  vortex  strength  and  size  or  shifting  of  the  vortex 
breakdown  location,  its  predicted  solutions  were  compared  qualitatively  rather  than 
quantitatively  with  experimental  data  from  UC  and  LSU. 
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Computational  Plan  of  Attack 

With  the  exception  of  the  sub-study  and  numerical  mesh  assessments  discussed 
above,  it  was  determined  that  blowing  cases  were  to  be  evaluated  using  a  time  accurate 
flow  solver  due  to  the  inherently  and  highly  unsteady  nature  of  the  vortical  flow  and 
breakdown.  All  blowing  cases  used  second-order  temporal  and  spatial  discretization  with 
CFL  =  5;  for  the  unsteady  cases,  At  =  0.0004  sec,  which  was  the  same  step-size  used  for 
non-blowing  cases.  Two  basic  strategies  were  attempted  with  the  unsteady  blowing 
cases:  initialize  and  run  the  case  with  blowing  active,  and  run  the  case  using  an  initial, 
fully-developed-flow  solution  from  a  no-blowing  case.  Both  strategies  proved 
unsuccessful  because  the  blowing  momenta  of  the  jets  were  numerically  overwhelmed  by 
the  primary  vortex  momentum  as  will  be  shown. 

Results 

Bottom  line  for  the  numerical  cases  with  along-core  blowing  is  that  all  three- 
dimensional  attempts  to  affect  vortex  breakdown  were  unsuccessful.  In  Figure  4.4, 
nondimensional  velocity  contours  from  a  two-dimensional,  steady  S-A  model  showed 
effects  from  blowing;  that  is,  vortex  breakdown  was  not  necessarily  eliminated,  but  its 
effects  were  clearly  changed  to  give  a  smaller  region  of  vortical  flow.  Interestingly,  the 
S-A  model  predicted  vortex  breakdown  in  the  two-dimensional  numerical  case  but  not  in 
the  three-dimensional  case.  Expanding  to  three  dimensions  introduced  numerical  issues 
which  were  not  resolved. 
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No  Blowing 


Port  1  Blowing  at  345  m/s 


Figure  4.4  -  Two-Dimensional  Solution  with  and  without  Blowing  at  Port  1  (Steady  S-A, 

a  =  15°,  Voo  =  15  m/s) 

Figure  4.5  shows  streamlines  (blue  lines)  and  vortex  core  centerline  identification 
(red  line)  from  Mesh  Revision  L’s  unsteady,  RSM  turbulence  cases  with  and  without 
blowing.  The  left  column  shows  the  predicted  solution  with  no  blowing  at  t  =  0.08,  0.16 
and  0.24  sec.  These  images  indicate  that  the  flow  solution  was  far  from  being  fully 
developed,  since  the  vortex  breakdown  position  continued  to  move  upstream  along  the 
wing  surface  over  time,  moving  from  x/c  =  0.75  to  about  0.55.  The  middle  column  of 
images  shows  the  predicted  solution  with  sonic  blowing  (345  m/s)  at  Port  1  -  all  other 
conditions,  models  and  solvers  being  the  same  as  those  from  the  left  column  -  and  the 
right  column  shows  the  results  with  sonic  blowing  at  all  three  ports.  The  case  with 
blowing  from  all  three  ports  was  attempted  in  an  effort  to  introduce  more  blowing 
momentum  into  the  primary  vortex  core.  Inspection  of  these  three  columns  of  images 
reveals  that  there  are  so  significant  differences  between  the  cases  at  the  same  time 
snapshots.  Each  predicts  vortex  breakdown  at  the  same  location,  indicating  failure  to 
delay  breakdown. 
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Figure  4.5  -  Vortex  Breakdown  via  Streamlines  for  Unsteady  RSM  Model  at  No  Blowing 
and  Blowing  Configurations  (Revision  L,  a  =  18°,  Vm=  16  m/s) 


Figure  4.6  shows  the  unsteady  RSM  solution  with  sonic  blowing  from  all  three 
ports  at  t  =  0.24  sec.  Nondimensional  velocity  contours  and  blue  streamlines  through  the 
vortex  core  indicate  that  vortex  breakdown  was  not  delayed  by  introducing  blowing 
momentum.  Nondimensional  velocity  vectors  were  included  to  show  that  each  port’s 
blowing  vector  was  precisely  aligned  with  the  vortex  core  centerline;  by  rotating  this 
graphic  and  inspecting  from  above  (image  not  included  here),  it  was  observed  that  the 
blowing  vector  lines  actually  bisected  the  core  streamlines.  Red  and  purple  streamlines 
emanate  from  the  blowing  ports  and  show  that  mass  flow  from  the  jets  was  swept  around 
the  vortex  instead  of  entering  its  core,  thus  negating  their  effect.  This  indicates  that  the 
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Figure  4.6  -  Effect  of  Port  1+2+3  Sonic  Blowing  on  Unsteady  RSM  Solution 
(Revision  L,  a  =  18°,  V00=  16  m/s) 


jet  momenta  were  either  overwhelmed  by  momentum  of  the  primary  vortex  or  dissipated 
before  reaching  the  vortex  core. 

Figure  4.7  compares  Cp  curves  at  the  four  chordwise  positions  for  cases  with  no 
blowing  and  with  sonic  blowing  from  Ports  1,  2  and  3,  using  Mesh  Revision  L,  unsteady 
RSM  turbulence  model  at  t  =  0.24  sec.  Blowing  resulted  in  no  notable  improvements  at 
x/c  =  0.35  and  0.95.  Peak  Cp  magnitudes  increased  by  8.2%  at  x/c  =  0.55  and  6.1%  x/c  = 
0.75.  Without  continuing  computation  of  this  blowing  case  to  roughly  t  =  1.6  sec,  it  is 
uncertain  whether  these  improvements  would  be  maintained.  Also  these  improvements 


4-10 


Baseline  vs  Blowing  at  345  m/s  from  Ports  1 ,2&3 ^ 
Unsteady  RSM  Turb  Model  with  Rev  L  at  t  =  0.24  sec 
left  to  right- x/c  =  0.35,  0.55, 0.75,  0.95 


y/s  (percent) 

l _ _ _ 

Figure  4.7  -  Cp  Comparison  of  No  Blowing  vs  Sonic  Blowing  at  Ports  1,  2,  and  3 
at  t  =  0.24  sec  (Revision  L,  Unsteady  RSM,  a  =  18°,  Vx  =  16  m/s) 


are  insignificant  in  comparison  with  those  obtained  in  experimental  testing  at  UC,  where 


peak  CP  magnitudes  increased  over  the  baseline  values  by  as  much  as  17,  1 1,  36  and  63% 


at  x/c  =  0.35,  0.55,  0.75  and  0.95,  respectively  (May,  2002. a). 


Figure  4.8  compares  contours  of  total  pressure  at  t  =  0.16  and  0.24  sec  for  cases 


with  no  blowing  and  with  sonic  blowing  from  Ports  1,  2  and  3,  using  the  unsteady  RSM 


turbulence  model.  These  images  show  that  the  vortex  core  remains  essentially  unaffected 


by  the  blowing.  The  images  from  the  blowing  case  (right  column)  show  “holes”  or  small 


regions  of  higher  pressure  downstream  of  the  blowing  ports  and  close  to  the  wing 


surface.  This  is  an  additional  indication  that  while  blowing  does  affect  the  flow,  it 


dissipates  prior  to  reaching  the  vortex  core  and  is  swept  into  the  vortex’  momentum. 
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No  blowing  Pori  1, 2  &  3  blowing  at  345  m/s 


Unsteady  RSM,  a  =  18°,  Vjnf  =  16  m/s 

Figure  4.8  -  Contours  of  Total  Pressure  for  Cases  with  No  Blowing  and  Sonic  Blowing  at 
All  Ports  (Revision  L,  Unsteady  RSM,  a  =  18°,  Vx  =  16  m/s) 

Since  altering  the  wing  angle  of  attack  would  likely  have  little  effect  on 

improving  this  issue  with  blowing  jet  momentum,  the  next  attempt  was  to  increase  the 

freestream  velocity  to  perhaps  ease  the  computational  complexity  of  predicting  flow  at 

low  freestream  velocity.  Figure  4.9  compares  results  from  four  cases  -  no  blowing  and 

blowing  with  Vx  =  16  m/s  (top  row);  and  no  blowing  and  blowing  with  V,,  =  103.8  m/s  or 

Moo  =  0.30  (bottom  row).  These  images  show  there  was  no  effect  from  Port  1  blowing  at 

two  different  freestream  velocities. 
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Figure  4.9  -  Effect  of  Port  1  Sonic  Blowing  on  Unsteady  RSM  Solution  for  M*  =  0.04 

and  0.3  (Revision  L,  a  =  18°) 


All  of  the  above  cases  initialized  blowing  at  t  =  0  sec.  In  another  sub-study  to 
show  any  effect  from  introducing  blowing  jet  momentum,  the  unsteady,  RSM  case  with 
blowing  from  Port  1  was  initialized  at  t  =  2.12  sec  from  the  case  with  no  blowing. 
Initializing  sonic  blowing  from  a  solution  with  fully  developed  flow  also  had  no  effect  on 
the  vortex  or  on  location  of  vortex  breakdown.  After  running  the  case  for  0.08  sec  (4,000 
iterations  or  200  time  steps),  the  solution  appeared  nearly  identically  to  the  top  right 
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image  in  Figure  4.9,  again  indicating  that  vortex  momentum  may  have  overwhelmed 


momentum  introduced  by  the  blowing  jet.  Predicted  Cp  and  Cl  values  changed  by  tenths 


of  a  percent,  where  the  LSU  study  showed  that  blowing  from  Port  1  increased  lift  by 


60%,  or  from  CL  =  0.8  to  1.28  (Guillot,  1999:  32,  50). 


In  another  case,  freestream  velocity  was  decreased  to  5  m/s,  where  blowing  was 
sonic  from  jet  Port  1,  using  a  steady  S-A  model.  An  attempt  with  zero  freestream 
velocity  was  numerically  unstable,  but  using  a  low  velocity  demonstrated  more  the 
inadequacy  of  the  jet  blowing  model  rather  than  the  possibility  of  numerically 
overwhelmed  flow  momentum.  As  with  Figures  4.6  and  4.9,  Figure  4.10  shows  that  the 
blowing  momentum  appeared  to  not  reach  the  vortex  structure  but  was  either 


immediately  swept  into  its  rotating  flow  or  quickly  dissipated  before  affecting  even  this 


tondim  Velocity  Contours  &  Vectors  (u_inf=5  m/s) 
slice  along  vortex  core  centerline 
Streamlines  thru  Core  and  from  Port  1 


Figure  4.10  -  Port  1  Sonic  Blowing  for  =  0.01  (Steady  S-A,  Revision  L,  a 
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low-energy  vortex.  With  this  low  freestream  velocity,  it  is  unlikely  that  the  jet  flow  was 
overwhelmed  by  the  vortical  structure;  this  is  perhaps  a  greater  indication  that  the  jet’s 
boundary  conditions  were  inadequately  specified  or  that  the  numerical  mesh  in  the  region 
of  blowing  was  not  sufficiently  refined.  Hence,  additional  investigation  into  the  jet 
blowing  profile,  boundary  condition  parameters,  blowing  region  mesh  refinement  to 
include  resolution  for  LES  and  secondary  vortex  is  merited. 

Thus  no  three-dimensional  attempts  in  this  numerical  study  successfully  predicted 
delay  of  vortex  breakdown  from  along-core  blowing  into  the  vortex. 
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V.  Conclusions 


Final  Configuration 

The  final  numerical  mesh  for  this  study  (Mesh  Revision  L)  modeled  the  half  delta 
wing  mounted  at  a  =  18°,  flush  against  a  boundary  layer  refresher  plate,  which  was  flush 
against  the  wind  tunnel  wall.  The  numerical  or  virtual  wind  tunnel  maintained  the  same 
dimensions  as  the  test  section  of  the  tunnel  used  during  testing  and  data  generation  at  the 
University  of  Cincinnati.  The  wing’s  angle  of  attack,  which  was  reported  to  be  15° 
during  physical  testing,  was  increased  in  the  numerical  model  to  give  better  correlation 
between  its  predicted  solution  and  experimental  data.  This  study  has  shown  that  there 
was  likely  error  in  the  reported  a.  Freestream  velocity  for  the  numerical  study  was  16 
m/s,  which  was  increased  within  an  acceptable  tolerance  of  the  reported  Vo0=  15.4  m/s  in 
an  effort  to  improve  data  correlation.  The  physical  boundary  layer  refreshing  plate  was 
offset  by  1 .27  cm  from  the  tunnel  wall,  but  this  study  showed  that  it  was  acceptable  to 
make  the  computationally  simplifying  assumption  of  placing  it  flush  against  the  tunnel 
wall. 

Numerical  results  were  investigated  and  compared  using  FLUENT  Version 
6.0.20’s  steady  and  unsteady,  second-order  accurate  flow  solvers,  and  the  following  flow 
models:  inviscid,  laminar,  and  turbulence  models,  Spalart-Allmaras  (S-A), 
Renormalization  Group  k-s,  Reynolds  Stress  (RSM),  and  Large  Eddy  Simulation  (LES). 
While  used  preliminarily,  the  LES  model  resulted  in  an  unacceptable  solution  because  the 
numerical  mesh  was  inadequately  refined  for  the  LES  wall  y+  requirement. 
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Conclusions 


This  study  revealed  no  combination  of  FLUENT  flow  solver  and  turbulence 
model  that  predicted  a  completely  acceptable  or  adequate  solution.  Time  averaged 
solutions  generally  converged  within  hours  and  predicted  quantitative  results  relatively 
close  to  experimental  data,  but  they  did  not  predict  most  of  the  unsteady  characteristics  of 
vortical  flow  and  vortex  breakdown.  Time  accurate  solutions  attained  fully  developed 
flow  typically  within  a  couple  weeks  and  only  the  RSM  turbulence  model  predicted  some 
of  the  unsteady  flow  physics;  further,  their  quantitative  predictions  were  not  an 
improvement  over  those  of  the  steady  solver.  Regarding  flow  models  in  unsteady  and/or 
steady  time,  laminar,  k-s,  RSM  and  LES  each  predicted  vortex  breakdown  but  not  to  an 
acceptable  quantitative  degree  for  vortex  size  and  strength;  the  S-A  turbulence  model 
gave  relatively  close  quantitative  prediction  of  the  solution  but  failed  to  predict  vortex 
breakdown.  A  recent  numerical  study  also  confirmed  that  the  S-A  model  failed  to 
predict  vortex  breakdown;  however,  that  study  also  revealed  that  modifications  to  the  S- 
A  model,  which  cater  to  flow  issues  associated  with  vortex  breakdown,  did  successfully 
predict  the  highly  unsteady  phenomenon  (Morton,  Forsythe,  Mitchell,  Hajek,  2002). 
Those  modified  S-A  models  are  not  currently  available  in  FLUENT. 

No  three-dimensional  numerical  case,  regardless  of  flow  solver,  turbulence  model 
or  freestream  velocity,  successfully  controlled  or  in  any  way  affected  vortex  breakdown. 
This  was  due  to  one  or  all  of  the  following:  along-core  jet  blowing  momentum  was 
numerically  overwhelmed  by  momentum  of  the  primary  vortex;  boundary  condition  for 
the  jet  was  inadequately  specified  and  resulted  in  rapid  dissipation  into  the  vortex 
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structure;  or  the  numerical  grid  in  the  regions  of  blowing  was  insufficiently  refined  to 
prevent  the  rapid  blowing  dissipation. 

The  numerical  grid  used  for  the  majority  of  cases  predicted  solutions  adequate  for 
comparative  analysis.  However,  the  mesh  was  not  sufficiently  refined  to  give  optimal 
solutions,  particularly  in  the  regions  over  the  forward  half  of  the  delta  wing,  along  the 
wing’s  leading  edge,  within  and  immediately  surrounding  the  primary  vortex,  and  in  the 
regions  of  jet  blowing. 

A  likely  candidate  for  most  of  the  numerical  problems  in  this  study  was  the  wing 
mounted  flush  against  the  tunnel  wall.  Since  there  were  no  data  for  comparison,  no 
numerical  study  was  performed  with  a  full  wing  model  inside  the  wind  tunnel.  That 
would  eliminate  impact  of  boundary  layer  flow  on  vortex  generation  and  breakdown  and 
is  the  current  focus  of  wind  tunnel  research  at  UC. 

Future  Work 

Using  a  numerical  mesh  or  developing  a  full  wing  numerical  model  with  greater 
resolution  in  the  regions  mentioned  above,  one  option  for  future  work  is  to  use  a  different 
flow  solver  software,  such  as  Cobalt.  Two  numerical  studies  of  a  delta  wing  at  high 
angle  of  attack  have  obtained  solutions  with  Cobalt  which  compared  closely  to 
experimental  data  and  which  adequately  predicted  the  unsteady  properties  of  vortex 
breakdown  (Cummings,  Morton,  Siegel,  2003;  Morton,  Forsythe,  Mitchell,  Hajek,  2002). 
Another  option  is  to  develop  three-dimensional  CFD  code  specific  to  this  problem  and 
which  implements  one  of  the  modified  S-A  transport  solvers  suggested  by  Morton  et  al. 
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If  better  computational  resources  become  available,  a  DNS  solver  would  be  best,  but  that 
is  likely  not  a  reasonable  option  for  many  years  hence. 

If  work  is  to  continue  with  the  FLUENT  flow  solver,  following  are  some 
recommendations.  Develop  a  full  delta  wing  computational  model,  including  wind 
tunnel  geometry  for  wall  effects;  UC  experimental  data  for  comparison  should  be 
available  within  the  year.  Refine  the  numerical  mesh  for  compatibility  with  the  unsteady 
LES  turbulence  model  and  divide  the  computational  load  among  a  sufficient  number  of 
processors  to  generate  a  solution  within  a  reasonable  amount  of  time.  Investigate 
FLUENT’S  other  turbulence  models  -  Realizable  k-s,  Standard  and  Shear-Stress 
Transport  k-co,  and  LES  with  RNG-Based  Sub-Grid  Scale  model  (FLUENT,  2001: 
10.4.3,  10.5.1-2,  10.7.2).  Dynamic  maneuvering  through  a  number  of  angles  of  attack 
may  be  simulated  by  creating  a  user-defined  function  for  the  pressure  farfield  boundary 
condition  (refer  to  FLUENT  6.0  UDF  Manual)',  this  would  require  modeling  the  wing  at 
a  =  0°  with  a  farfield  boundary.  As  part  of  troubleshooting  the  blowing  cases,  create  a 
user-defined  function  for  the  velocity  profile  at  the  blowing  port  boundaries,  giving  it  the 
more  correct  parabolic  form;  consider  significantly  greater  refinement  of  the  numerical 
mesh  in  the  blowing  regions;  and  conduct  more  research  into  blowing  jet  boundary 
conditions. 
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Appendix  A:  Experimental  Data 


Following  are  the  tabulated  experimental  data  obtained  during  wind  tunnel  testing  at 
University  of  Cincinnati  (UC).  Ambient  conditions  were  reported  at  T  =  25°  C  (298  K) 
and  p  =  29.4  in.  Hg  (99.56  kPa).  Wing  a  =  15°  and  F«  =  15.4  m/s.  Baseline  data, 
without  along-core  blowing,  were  sampled  at  100  Hz  for  a  period  of  5  sec  then  averaged 
for  the  values  in  Table  A.l.  As  no  numerical  model  predicted  the  effects  of  along-core 
blowing  into  the  vortex  core,  those  corresponding  experimental  data  are  not  provided. 


Table  A.l 


Baseline  Cp  Data  over  Wing  Surface,  N 


o  Blowing 


x/c 

y/s 

sensor 

static  P  (Pa) 

dyn  P  (Pa) 

CP 

0.35 

0.15 

1 

-58.516 

139.683 

0.419 

0.17 

2 

-67.633 

139.993 

0.483 

0.19 

3 

-97.310 

139.886 

0.696 

0.20 

4 

-166.452 

139.886 

1.190 

0.22 

5 

-273.960 

140.121 

1.955 

0.25 

6 

-314.272 

140.110 

2.243 

0.55 

0.25 

7 

-43.956 

140.228 

0.313 

0.29 

8 

-61.847 

140.345 

0.441 

0.31 

9 

-99.160 

140.665 

0.705 

0.33 

10 

-169.803 

140.313 

1.210 

0.36 

11 

-232.605 

140.750 

1.653 

0.39 

12 

-229.654 

141.775 

1.620 

0.41 

13 

-212.346 

140.857 

1.508 

0.43 

14 

-195.944 

141.113 

1.389 

0.75 

0.35 

15 

-34.196 

141.487 

0.242 

0.41 

17 

-86.265 

141.081 

0.611 

0.45 

18 

-133.785 

141.839 

0.943 

0.48 

19 

-173.874 

141.668 

1.227 

0.51 

20 

-181.919 

141.263 

1.288 

0.55 

21 

-162.829 

141.850 

1.148 

0.58 

22 

-139.113 

141.199 

0.985 

0.95 

0.41 

23 

-37.118 

141.540 

0.262 

0.45 

24 

-46.050 

142.074 

0.324 

0.50 

25 

-64.438 

141.807 

0.454 

0.54 

26 

-83.499 

141.529 

0.590 

0.58 

27 

-100.232 

142.138 

0.705 

0.62 

28 

-108.530 

142.010 

0.764 

0.66 

29 

-107.235 

141.924 

0.756 

0.70 

30 

-93.804 

142.362 

0.659 
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Since  the  dynamic  pressure  appeared  to  fluctuate  in  these  data,  the  CP  values 
were  normalized  for  comparison  with  computational  data  by  providing  a  constant 
dynamic  pressure,  based  on  lab  ambient  density  and  tunnel  freestream  velocity  -  such 
that  dynamic  pressure  became  143.206  Pa.  Table  A. 2  shows  the  data  which  were  then 
used  to  compare  with  numerically  predicted  values. 


Table  A.2  -  Normalized  Base' 


ine  Cp  Data  over  Wing  Surface,  No  Blowing 


x/c 

y/s 

sensor 

static  P  (Pa) 

dyn  P  (Pa) 

CP 

0.35 

0.15 

1 

-58.516 

143.206 

0.409 

0.17 

2 

-67.633 

143.206 

0.472 

0.19 

3 

-97.310 

143.206 

0.680 

0.20 

4 

-166.452 

143.206 

1.162 

0.22 

5 

-273.960 

143.206 

1.913 

0.25 

6 

-314.272 

143.206 

2.195 

0.55 

0.25 

7 

-43.956 

143.206 

0.307 

0.29 

8 

-61.847 

143.206 

0.432 

0.31 

9 

-99.160 

143.206 

0.692 

0.33 

10 

-169.803 

143.206 

1.186 

0.36 

11 

-232.605 

143.206 

1.624 

0.39 

12 

-229.654 

143.206 

1.604 

0.41 

13 

-212.346 

143.206 

1.483 

0.43 

14 

-195.944 

143.206 

1.368 

0.75 

0.35 

15 

-34.196 

143.206 

0.239 

0.41 

17 

-86.265 

143.206 

0.602 

0.45 

18 

-133.785 

143.206 

0.934 

0.48 

19 

-173.874 

143.206 

1.214 

0.51 

20 

-181.919 

143.206 

1.270 

0.55 

21 

-162.829 

143.206 

1.137 

0.58 

22 

-139.113 

143.206 

0.971 

0.95 

0.41 

23 

-37.118 

143.206 

0.259 

0.45 

24 

-46.050 

143.206 

0.322 

0.50 

25 

-64.438 

143.206 

0.450 

0.54 

26 

-83.499 

143.206 

0.583 

0.58 

27 

-100.232 

143.206 

0.700 

0.62 

28 

-108.530 

143.206 

0.758 

0.66 

29 

-107.235 

143.206 

0.749 

0.70 

30 

-93.804 

143.206 

0.655 
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Appendix  B:  Traub’s  Simple  Prediction  of  Vortex  Breakdown  Location 


Following  is  Traub’s  sequence  of  equations  used  to  predict  location  of  vortex 
breakdown.  Since  this  model  is  based  upon  curve-fitting  to  experimental  data  for  delta 
wings  with  sweep  angle,  A  =  65,  70,  75  and  80°,  the  results  in  this  case  for  A  =  60°  were 
only  valid  for  ensuring  the  predicted  results  were  within  an  acceptable  range  or 
“ballpark”  estimate  of  what  they  should  be  (Traub,  1996). 

The  first  of  three  steps  was  to  determine  aBD_TE ,  the  angle  of  attack  for  a  delta 
wing  at  which  vortex  breakdown  occurs  at  the  trailing  edge: 

tan  aBD_TE  =13.47  •  tans  ■  e~6'9'£  (B.l) 

where  a  is  wing  apex  half  angle  or  0.5  A.  For  A  =  60°,  aBD_TE  =  1 1 .848°. 

Second,  the  nondimensional  circulation  must  be  determined.  It  is  defined  as 

r 

- =  4.63  •  tan0  s  s  •  tan1'2  a  •  cosa  (b.2) 

^00  '  Cf 

where  r  is  circulation,  a  is  wing  angle  of  attack,  cr  is  wing  root  chord  length,  and  Vx  is 
freestream  velocity.  This  nondimensional  quantity  must  be  determined  using  a  =  aBD  TE  , 
as  well  as  using  the  angle  of  attack  of  interest.  For  the  ocBD_TE  determined  above  and  for 
a  =  15  and  18°,  the  nondimensional  circulation  values  become  0.44824,  0.59339  and 
0.73633,  respectively. 

Finally,  the  predicted  breakdown  location  is  detennined  by 


B-l 


(B-3) 


r 

1  BD-TE 

r 

A  a>BD-TE 

where  n  is  an  empirically  determined  curve-fit  constant.  Traub  used  n  =  3  to  fit  the  data 
for  A  =  65°,  so  that  value  was  also  used  for  this  case  with  A  =  60°.  Thus  the  final  answer 
or  predicted  locations  for  vortex  breakdown  for  A  =  60°  at  a  =  15  and  18°  were  0.431  and 
0.226,  respectively. 
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Appendix  C:  Gridgen  Lessons 


Following  are  some  lessons  learned  from  creating  a  numerical  mesh  using  the 
Gridgen  software  and  from  interfacing  Gridgen  and  FLUENT.  One,  the  user  must  follow 
the  default  directional  indicators  in  the  graphical  user  interface  display  when  selecting 
lines  to  create  borders  for  a  domain,  else  it  may  result  in  a  negative  volume  for  the 
associated  block;  FLUENT  will  not  accept  a  mesh  with  negative  volumes.  An  exception 
to  this  is  for  creating  a  domain  with  an  embedded  form,  where  the  user  must  select  lines 
in  the  direction  opposite  to  that  used  for  the  domain  exterior.  Two,  when  declaring 
boundary  conditions  for  a  problem  with  confined  flow,  if  any  of  them  are  inflow  (e.g., 
pressure  gradient,  velocity,  or  mass  flow  rate,  and  not  farfield  pressure),  the  exit  must  be 
an  outlet  (pressure,  velocity,  or  mass  flow)  rather  than  a  simple  outflow  condition.  Three, 
when  interfacing  structured  and  unstructured  volumes,  Gridgen  creates  “Type  I” 
boundaries  by  default;  the  user  should  not  specify/change  that  boundary  condition  to  be 
an  interface,  else  FLUENT  will  not  correctly  import  it. 

Last,  if  two  or  more  boundaries  within  a  given  block  have  the  same  boundary 
condition,  FLUENT  will  automatically  group  them,  such  that  they  must  then  have  the 
same  initial  condition.  For  example,  if  two  blowing  jets  on  the  wing  surface  are  specified 
as  velocity  inlet  boundary  conditions,  they  must  both  be  active  or  inactive,  both  have  the 
same  blowing  angle  and  velocity,  and  so  on.  Similarly  for  another  example,  if  the  wing 
and  tunnel  walls  are  all  specified  with  a  wall  boundary  condition,  FLUENT  will  group 
them  and  compute  quantities  such  as  CL  and  CD  based  on  the  entire  surface  area  and  not 
just  the  wing.  To  avoid  this  pitfall,  in  Gridgen  the  user  must  specify  related  but  different 
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boundary  conditions  for  any  boundary  that  needs  to  be  distinct;  then  in  FLUENT,  the 
user  may  specify/change  those  boundaries  to  what  they  need  to  be. 
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Appendix  D:  Matrix  of  Transformation  to  Determine  Jet  Blowing  Angles 


Following  is  the  development  of  the  matrix  of  transformation  used  to  detennine 
correct  blowing  angles  for  the  momentum  jets  on  the  half  delta  wing’s  upper  surface. 
This  transfonnation  allows  for  any  variation  in  orientation  of  the  wing  or  its  blowing  jets. 
The  matrix  in  Equation  D.l  uses  Euler  angles  for  yaw  ( if/),  pitch  ( 6 )  and  roll  ((f))  to 
transform  body  (wing)  reference  frame  to  a  Cartesian  fixed  reference  frame  (in  that  order 
of  rotation),  which  is  then  used  as  FLUENT  boundary  conditions  for  the  blowing  jets. 


• 

Xf 

u 

• 

yf 

• 

=[4 

V 

Zf 

w 

(D.l) 


where  Xf-dot,  Vf-dot  and  Zf-dot  are  fixed  reference  Cartesian  velocity  components,  A  is 
defined  as 


cosO-cos^  sin^-sin#-cos^-cos^-sin^  cos^-sin#-cos^  +  sin^-sin^ 
cos#- sin  ^  sin^-sin#-sin^  +  cos^-cos^  cos^-sin#-sin^-sin^-cos^ 
-sin#  sin^-cos#  cos^-cos# 


(D.2) 


and  u,  v  and  w  are  three-dimensional  components  of  the  blowing  velocity  defined  by 

u  =  -Vjet-cos6ercosOaz  (D.3) 


v  ~  V jet '  cos#e/  •  sin#az 


w=-Vjet-sm0ei 


(D.4) 

(D.5) 
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where  Vjet  is  velocity  magnitude  of  blowing  jet,  9ei  is  elevation  angle  of  blowing  jet,  and 
Qaz  is  azimuthal  angle  of  blowing  jet  (measured  counter-clockwise  from  a  line  parallel  to 
the  wing’s  root,  as  shown  in  Figure  2.2);  the  blowing  jets  were  physically  limited  to 
variation  in  these  two  degrees  of  freedom.  Equation  D.2  came  from  Nelson  (1998:  102). 
Combining  Equations  D.l-5  gives  the  necessary  Cartesian  jet  blowing  components. 

For  this  case,  (f)  =  y/  =  0°,  and  pitch  angle  9  was  negative,  due  to  orientation  of 
the  wing’s  coordinate  axes.  Thus  Equations  D.l-5  simplify  to 

-  cos  cost?,,  -cos(?0_  -sinO-sint?,, 

=  viet  ■  cos0e/-sin0az  (D.6) 

sin 0  •  cos 6e]  •  cos Gaz  -  cos^  •  sin^ 
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